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FOREWORD 


This report is a revision of the document “MASTERFIT — 1987”, dated December 15, 1987, 
which it supersedes. A number of model revisions and improvements were made during 1988-91. 
They are briefly enumerated in the abstract. The computer code was also considerably revised during 
1988-91 to facilitate solution of large-scale problems. The new software still adheres to the basic 
MASTERFIT structure but, to prevent confusion concerning practical details, is named MODEST 
(for MODel and ESTimate). The present document corresponds to MODEST version 137, which has 
been in use since June, 1991. The author hopes to publish revisions of this document in the future, 
as modeling improvements warrant. 
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ABSTRACT 


This report is a revision of the document “MASTERFIT - 1987”, dated December 15, 1987, 
which it supersedes. Changes during 1988-91 included introduction of the octupole component of 
solid Earth tides, the NTJVEL tectonic motion model, partial derivatives for the precession constant 
and source position rates, the option to correct for source structure, a refined model for antenna 
offsets, modeling the unique antenna at Richmond, Florida, improved nutation series due to Zhu, 
Groten, and Reigber, and reintroduction of the old (Woolard) nutation series for simulation purposes. 
Text describing the relativistic transformations and gravitational contributions to the delay model has 
also been revised in order to reflect the computer code more faithfully. 
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SECTION 1 


INTRODUCTION 

In applications of radio interferometry to geodynamics and astrometry, observed values of delay 
and delay rate obtained from observations of many different radio sources must be passed simulta- 
neously through a multiparameter estimation routine to extract the significant model parameters. As 
the accuracy of radio interferometry has improved, increasingly complete models for the delay and 
delay rate observables have been developed. This report describes the current status of the delay 
model used in the Jet Propulsion Laboratory multiparameter estimation program “MODEST”, which 
is the successor to the “MASTERFIT” code developed at JPL in the 1970s. It is assumed that the 
reader has at least a cursory knowledge of the principles of VLBI. Some references which provide an 
introduction are the book by Thompson, Moran, and Swenson (1986), and two reports by Thomas 
(1981, 1987). 

The delay model is the sum of four major model components: geometry, clock, troposphere, 
and ionosphere. Sections 2 through 5 present our current models for these components, as well as 
their partial derivatives with respect to parameters that are to be adjusted by multiparameter fits 
to the data. The longest section (2) deals with the purely geometric portion of the delay and covers 
the topics of time definitions, tidal and source structure effects, coordinate frames, Earth orientation 
(universal time and polar motion), nutation, precession, Earth orbital motion, wave front curvature, 
gravitational bending, and antenna offsets. Section 6 describes the technique used to obtain the delay 
rate model from the delay model. Section 7 gives the values of physical constants used in MODEST, 
while section 8 outlines model improvements that may be required by more accurate data in the 
future. 
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SECTION 2 


GEOMETRIC DELAY 

The geometric delay is that interferometer delay which would be measured by perfect instrumen- 
tation, perfectly synchronized, if there were a perfect vacuum between the observed extragalactic or 
Solar-System sources and the Earth-based instrumentation. For Earth-fixed baselines, this delay can 
be as large as 20 milliseconds, changing rapidly (by up to 1.5 /i sec per second) as the Earth rotates. 
In general the geometric component is by fax the largest component of the observed delay. The main 
complexity of this portion of the model arises from the numerous coordinate transformations necessary 
to relate the reference frame used for locating the radio sources to the Earth-fixed reference frame in 
which station locations are represented. 

In the following we will assume, unless otherwise stated, that “celestial reference frame” means 
a reference frame in which there is no net proper motion of the extragalactic radio objects which 
are observed by the interferometer. This is only an approximation to some truly “inertial” frame. 
Currently, this celestial frame implies a geocentric, equatorial frame with the equator and equinox of 
J2000 as defined by the 1976 IAU conventions, including the 1980 nutation series (Seidelmann, 1982, 
and Kaplan, 1981). 

In this equatorial frame, some definition of the origin of right ascension must be made. We will 
not discuss that in this report, since one definition is at most a rotation from some other definition, 
and can be applied at any time. The important point is that consistent definitions must be used 
throughout the model development. The need for this consistency will, in all probability, eventually 
lead to our defining the origin of right ascension by means of the JPL planetary ephemerides, followed 
by our using interferometric observations of both natural radio sources and spacecraft at planetary 
encounters as a means of connecting the planetary and the radio reference frames (Dewey, 1991, 
Newhall et al., 1986). 

Also, unless otherwise stated, we will mean by “terrestrial reference frame” some reference frame 
tied to the mean surface features of the Earth. Currently, we are using a right-handed version of the 
CIO reference system with the pole defined by the 1903,0 pole. In practice, this is accomplished by 
defining the position of one of the interferometric observing stations (generally DSS 14 at the Goldstone 
Deep Space tracking complex), and then by measuring the positions of the other stations under a 
constraint. This constraint is that the determinations of Earth orientation agree on the average with 
the International Earth Rotation Service (IERS) (1991) [and its predecessor, Bureau International 
de 1’Heure (BIH) (1983)] measurements of the Earth’s orientation over some substantial time interval 
(» years). This procedure, or its functional equivalent, is necessary since the interferometer is sensitive 
only to the baseline vector as measured in the celestial frame. The VLBI technique does not have 
any preferred origin relative to the structure of the Earth. The rotation of the Earth does, however, 
provide a preferred direction in space which can be associated indirectly with the surface features of 
the Earth. 

In contrast, geodetic techniques which involve the use of artificial satellites, or the Moon, are 
sensitive to the center of mass of the Earth as well as the spin axis. Thus, those techniques require 
only a definition of the origin of longitude. We anticipate that laser ranging to the retroreflectors on 
the Moon (LLR) will allow a realizable practical definition of a terrestrial frame, accurately positioned 
relative to a celestial frame which is tied to the planetary ephemerides. The required collocation of 
the laser and VLBI stations is being provided by Global Positioning Satellite (GPS) measurements 
of baselines between VLBI and laser sites starting in the late 1980s (e.g., Ray et al 1991). Careful 
definitions and experiments of this sort will be required to realize a coordinate system of centimeter 
accuracy. In the meantime, we must establish interim coordinate systems carefully enough so that we 
do not degrade the intrinsic accuracy of the interferometer data by introducing “model noise”. 

The relativistic delay formulation presented in this report is the same as that in an earlier report 
(Sovers and Fanselow, 1987) except for a small change in the gravitational correction. Among the 
estimated parameters, only baseline length is affected by this change, in that all distances are increased 
by the same factor of « 2 parts in 10®. Special relativistic terms in the model delay have not been 
changed from the earlier report. 


2 



Except for subcentimeter relativistic complications caused by the locally varying Earth potential 
(as discussed below), calculation of the VLBI model for the observed delay can be summarized as 
follows: 

1 Specify the proper locations of the two stations as measured in an Earth-fixed frame at the time 
that the wave front intersects station #1. Let this time be the proper time t\ as measured by a 

clock in the Earth-fixed frame. _ . 

2. Modify the station locations for Earth-fixed effects such as solid Earth tides, tectonic motion, 

and other local station motion. _ , , 

3. Transform these proper station locations to a celestial coordinate system with the origin at the 
center of the Earth, but moving with the Earth. This is a composite of 10 separate rotations, 

represented by a rotation matrix Q(t ). . , 

4. Perform a Lorentz transformation of these proper station locations from the geocentric celestial 
frame to a frame at rest relative to the center of mass of the Solar System, and rotationally 
aligned with the celestial geocentric frame. 

5. In this Solar-System-barycentric frame, compute the proper time delay for the passage of the 
specified wave front from station #1 to station #2. Correct for source structure. Also, add in the 
effective change in proper delay caused by the differential gravitational retardation of the signal. 

6. Perform a Lorentz transformation of this SSB geometric delay back to the celestial geocentric 
frame moving with the Earth. This produces the adopted model for the geometric portion of the 
observed delay. 

7. To this geometric delay, add the contributions due to clock offsets, to tropospheric delays, and 
to the effects of the ionosphere on the signal (see sections 3 through 5). 

As indicated in step 5, the initial calculation of delay is carried out in a frame at rest relative 
to the center of mass of the Solar System (SSB frame.) First, however, steps 1 through 4 are carried 
out in order to relate proper locations in the Earth-fixed frame to corresponding proper locations 
in the SSB frame. Step 4 in this process Lorentz transforms station locations from the geocentric 
celestial frame to the SSB frame. This step incorporates special-relativistic effects to all orders of v/c. 
In the presence of gravity, this transformation can be viewed as a special relativistic transformation 
between proper coordinates of two local frames (geocentric and SSB) in relative motion. For both 
frames, the underlying gravitational potential can be viewed approximately as the sum of locally 
constant potentials caused by all masses in the Solar System. The complications caused by small local 
variations in the Earth’s potential are discussed below. Initial proper delay is then computed (step 5) 
in the SSB frame on the basis of these SSB station locations and an a priori SSB source location. A 
small proper-delay correction is then applied to account for the differential gravitational retardation 
introduced along the two ray paths through the Solar System, including retardation by the Earth s 
gravity. A final Lorentz transformation including all orders of v/ c then transforms the corrected SSB 
proper delay to a model for the observed delay. 

Since the Earth’s potential varies slightly across the Earth (AC/e/c 2 « 4 X 10 10 from center 
to surface), the specification of proper distance is not as straightforward with respect to the Earth s 
potential as it is with respect to the essentially constant potentials of distant masses. To overcome this 
difficulty, output station locations are specified in terms of the “TDT spatial coordinates” (Shahid- 
Saless et al ., 1991) used in Earth-orbiter models. Baselines modeled on the basis of this convention 
deviate slightly in length (< 2 cm) from the proper values. A proper length that corresponds to a 
modeled baseline can be obtained through appropriate integration of the local metric (Shahid-Saless 
et al, 1991). In practice, such a conversion is not necessary since comparison of baseline measurements 
obtained by different groups would be carried out in terms of TDT spatial coordinates. 

The current model has been compared (Thomas, 1991, Treuhaft, 1991) with the a l-picosecond 
relativistic model for VLBI delays developed by Shahid-Saless et al. (1991). When reduced to the same 
form, the model presented here is identical to that model at the picosecond level, term by term, with 
one exception. Treuhaft and Thomas (1991) show that a correction is needed to the Shahid-Saless 
et al SSB system modeling of the atmospheric delay. This correction changes the Shahid-Saless et al 
result by as much as 10 picoseconds. The remainder of this section provides the details for the first 
six steps of the general outline above. 
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2.1 TIME INTERVAL FOR THE PASSAGE OF A WAVE FRONT 
BETWEEN TWO STATIONS 

The fundamental part of the geometric model is the calculation (step # 5 above) of the time 
interval for the passage of a wave front from station #1 to station # 2 . We actually do that calculation 
in a coordinate frame at rest relative to the center of mass of the Solar System. This part of the model 
is presented first to provide a context for the subsequent sections, all of which are heavily involved 
with the details of time definitions and coordinate transformations. We will use the same subscript 
and superscript notation which is used in section 2.7 to refer to the station locations as seen by an 
observer at rest relative to the center of mass of the Solar System. 

First, we calculate the proper time delay that would be observed if the wave front were planar. 
Next, we generalize this calculation to a curved wave front, and finally, we take into account the 
incremental effect which results from the fact that we must consider wave fronts that propagate 
through the various gravitational potential wells in the Solar System. 

2.1.1 Plane Wave Front 



IT AT TIME t. 

Figure 1. Geometry for calculating the transit time of a plane wave front 


Consider the case of a plane wave moving in the direction, k, with station 2 having a mean 
velocity, ^8 2 > as shown in figure 1. As mentioned above, distance and time are to be represented as 
proper coordinates in the SSB frame. The speed of light, which is c in this representation, is set equal 
to 1 in the following formulation. The proper time delay is the time it takes the wave front to move 
the distance l at speed c . This distance is the sum of the two solid lines perpendicular to the wave 
front in figure 1: 

t* 2 ~ ti = k ■ [r 2 («i) -r x (t!)] + k ^2 ~ f i] (2.1) 

This leads to the following expression for the geometric delay: 


_ k ■ [r 2 (t x ) 

1-k /9 2 
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The baseline vector, r 2 (ti) is computed on the basis of proper station locations calculated 

according to Eq. (2.155) below. 

2.1.2 Curved Wave Front 

In the case of a signal generated by a radio source within the Solar System it is necessary to 
include the effect of the curvature of the wave front. As depicted in figure 2, let a source irradiate two 
Earth-fixed stations whose positions are given by r*(t) relative to the Earth’s center. The position of 
the Earth’s center, R c (*i), as a function of signal reception time, t 1} at station #1 is measured relative 
to the position of the emitter at the time, t e , of emission of the signal received at time t x . While this 
calculation is actually done in the Solar System barycentric coordinate system, the development that 
follows is by no means restricted in applicability to that frame. 



Suppose that a wave front emitted by the source at time t c reaches station #1 at time and 
arrives at station #2 at time The geometric delay in this frame will be given by: 

T = t* 2 -t 1 = \R 2 (t* 2 )\-\R 1 (t 1 )\ (2.3) 

where all distances are again measured in units of light travel time. If we approximate the velocity of 
station #2 by 

92 - tl 

and use the relation 

Ri(ti) = R c (tj) + i\(ii) 


(2.4) 

(2.5) 
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we obtain: 

r = |R c (fi) 4 r 2 (ti) 4- P 2 A — iRc(ti) 4 r i(*i)l 



— Rc(h) [ |R-c + * 2 ] — |®-c + *x| ] 

(2.6) 

where 

r 2 (<l) +P 2 T 

* 2 Rc{t l) 

(2.7) 

and 

r l(«l) 

1 Rc(t l) 

(2.8) 

For e\ and e 2 < 10" 4 , we need to keep only terms of order e 3 in a sixteen-place machine in 
expand the expression for r in equation (2.6). This gives us: 

order to 


R c [** 2 (^ 1 ) r i(^i)] j i2 c A c (r) 

[l-Rc/M 2[1-R c /? 2 ] 

(2.9) 

where to order e 3 



A c(r) = [el-e\]- 

[(R c « 2 ) 2 + (R c *i) 2 + (Rc * 2 ) 3 - (Rc *2)e 2 - (R c *i) 3 + (Rc *i) e i] 

(2.10) 


The first term in (2.9) is just the plane wave approximation, t.e., as i? c — > oo, R c — ♦ k, with the 
second term in brackets in (2.10) approaching zero as r 2 /R c . Given that the ratio of the first term 
to the second term is w r/R ct wave front curvature is not calculable in a sixteen-place machine for 
R > 10 16 X r. For Earth-fixed baselines that are as long as an Earth diameter, requiring that the 
effects of curvature be less than 0.01 cm implies that the above formulation (2.10) must be used for 
R < 1.4 X 10 15 km, or approximately 150 light years. 

The procedure for the solution of (2.9) is iterative for e < 10 -4 , using the following: 


T n 


T 0 4 


Rc^cj^n— l) 

2[1 - R*£ 2 ) 


where 


Tq — Tpi ane wave 


For e > 10“ 4 , directly iterate on the equation (2.6) itself, using the procedure: 


( 2 . 11 ) 


( 2 . 12 ) 


r n — i? c |R c 4 * 2 ( r »-i) | ~ # c |Rc 4 *i| 


(2.13) 


where again tq is the plane wave approximation. 


2.1.3 Gravitational Delay 


Because a light signal propagating in a gravitational potential is retarded relative to its motion 
in field-free space, the computed value for the differential time of arrival of the signals at ri(ti) and 
r 2 (t 2 ) must te corrected for gravitational effects. For the geometry illustrated in figure 3, the required 
correction to coordinate time delay is given by Moyer (1971) as: 


A (l 4 'Tppjv )a*p 

— 3 


r $ 4 r 2 (t 2 ) 4 r, 2 
r» 4 r 2 (t 2 ) - t 92 


r 9 4 ri(tx) 4 r 9l 
r B 4 ri(ti) - r,i 


where r 9 i is defined as: 


t 9% — | r t(^t) r a (t c )| 


(2.14) 


(2.15) 
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Here i PrN 
1973): 


is the 7 factor in the parametrized post-Newtonian gravitational theory (e.g. Misner et al., 


1 + w 
1ffn = 2 + w 


(2.16) 


where w is the coupling constant of the scalar field. For general relativity, 7 p PJV = 1> w ~ 
However, we allow ^ PPN to be an estimated parameter so that by setting i PPN = -1, we also have 
the option of “turning off” the effects of general relativity on the estimate of the delay. This proves 
useful for software development. The gravitational constant, fi p} is 


Mp = Gm p 


(2.17) 


where G is the universal gravitational constant, and m p is the mass of the pth gravitating body. 



Figure 3. A schematic representation of the geodesic connecting two points in the 
presence of a gravitational mass 


Dropping the time arguments in (2.14), we have: 


Ag p = 


(l~f~ 7pfN )Mp 


r, + r 2 + r, 2 

r, + ri - r,i 

r 9 + ri + r 9 \ 

r, + r 2 - r , 2 J 


(2.18) 


This formulation is fine for r, « r< w r sii but can be put in a computationally better form for the 
case of distant sources with closely spaced VLBI receivers, i.e., |r 2 — r i|/ri — ► 0 , r*/r, — ► 0 . For these 
sources, expand Ag p in terms of Vijr 9} r 9 i/r 9i an< ^ make use of the relationship 


2r. 


■Ti + r, 2 ] 1/2 « r. 


Ti r. 


(2.19) 


This leads to 


A d - ppn)^p 

A Gp= ^ 


In 


ri +r r r, 


r 2 + r 2 • r. 


( 2 . 20 ) 


7 



( 2 . 21 ) 


for r x /r, — *■ 0 . 

If we further require that |r 2 — r x |/r x — ► 0 , and make use of 


r 2 = r x + Ar 


then: 


r 2 


r *21 

+ r 2 ■ r, = r x |^1 + 2 r i • A r/r x + (Ar/r x ) J + r x - r, + Ar • r, 
« r x (l + f x • Ar /r x ) + r x r, + Ar f. 

In the limit of Ar/ri — ► 0 : 

r 2 (1 +? 2 *?«) ri(l+fi r,) 4- Ar • (?i 4-r,) 
Substituting into (2.20) and expanding the logarithm, we obtain: 

(l + 7 r ™K III ~ r i) ' (?i + ? ») 


A Gp = -- 


r x (l + r x •?,) 


( 2 . 22 ) 


(2.23) 


(2.24) 


Using whichever of these three formulations (2.18, 2.20 or 2.24) is computationally appropriate, 
the model calculates a correction A c p for each of the major bodies in the Solar System (Sun, planets, 
Earth, and Moon). 

Before the correction A Gp can be applied to a proper delay computed according to Eq. ( 2 . 2 ), it 
must be converted from a coordinate-delay correction to a proper-delay correction appropriate to a 
near-Earth frame. For such proper delays, the gravitational correction is given to good approximation 
by 

A'g p = A Gp -(1 + 1 „s)Ut (2.25) 

where r is the proper delay given by Eq. ( 2 . 2 ), and where U is the negative of the gravitational 
potential of the given mass divided by c 2 , as observed in the vicinity of the Earth (U is a positive 
quantity). The Ur term is a consequence of the relationship of coordinate time to proper time, and 
the 7 ppn Ut term is a consequence of the relationship of coordinate distance to proper distance. 

The total gravitational correction used is: 


N 

&G = X/ ^ g p (2.26) 

p~ i 

where the summation to N is over the major bodies in the Solar System. For the Earth, the 
(1 4- 7 ppn )Ut term in Eq. (2.25) is omitted if one wishes to conform with the “TDT spatial coordi- 
nates” used to reduce Earth-orbiter data. The scale factor ( 14*7 pp N )U is approximately 1.97 x 10 -8 
for the Sun. A number of other conventions are possible. One of these, which does not omit the 
(1 4- 7 ppn)Ut term for the Earth, but evaluates it at the Earth’s surface, yields an additional scale 
factor of 0.14 X 10 -8 . In either case, the model delay is decreased. Consequently, all inferred “mea- 
sured” lengths increase by the same fraction relative to previous lengths ( t.g . by 19.7 parts per billion 
or 21.1 ppb). 

Some care must be taken in defining the positions given by r a , an< 3 ri(ti). We have chosen 

as the origin the position of the gravitational mass at the time of closest approach of the received 
signal to that object. The position, r 81 of the source relative to this origin is the position of that 
source at the time, t e , of the emission of the received signal. Likewise, the position, r,(t*), of the tth 
receiver is its position in this coordinate system at the time of reception of the signal. Even with 
this care in the definition of the relative positions, we are making an approximation, and implicitly 
assuming that such an approximation is no worse than the approximations used by Moyer (1971) to 
obtain (2.14). 
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Some considerations follow, regarding the use of appropriate times to obtain the positions of the 
emitter, the gravitational object, and the receivers. For a grazing ray emitted by a source at infinity, 
using the position of the gravitating body G at the time of reception of the signal at station ^1 rather 
than at the time of closest approach of the signal to G can cause a 15-cm error on baselines with 
a length of one Earth radius as shown by the following calculation. From figure 4, the calculated 
distance of closest approach, R } changes during the light transit time, tn g ht transit , of a signal from a 
gravitational object at a distance Reg by: 

A R « R EG e ■ t light tran , H = 0 • Rla/c (2.27) 


EMITTER 



Figure 4. A schematic representation of the motion of a gravitating object during the 
transit time of a signal from the point of closest approach to reception by 
an antenna 


Since the deflection is: 


A© « 2 


(l ~h ~ippN )a*p 


S{ A©) - -A© 


rARi 
. R - 


- A© 


cReg © 



Reg 1 

c© J dt 


(2.28) 

(2.29) 
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We consider the two bodies of largest mass in the Solar System: the Sun and Jupiter. For grazing 
rays, their respective deflections A0 are 8480 and 73 nanoradians. The barycentric angular velocities 
^ are estimated to be 0.06 and 17 nrad/sec for the Sun and Jupiter. Note that Eq. (2.27) does not 
apply to the Sun. The Sun’s motion in the barycentric frame has a period of 11 years with a radius 
of the order of the Sun’s radius. Using approximate radii and distances from Earth to estimate Reg 
and 0, Eq. (2.29) gives 25 nrad for Jupiter; the corresponding value for the Sun is 0.07 nrad. For 
a baseline whose length equals the radius of the Earth, 6(AS)R E is thus approximately 0.05 and 15 
cm for the Sun and Jupiter, respectively. The effect is much smaller for the Sun in spite of its much 
larger mass, due to its extremely slow motion in the barycentric frame. 

In view of the rapid decrease of gravitational deflection with increasing distance of closest ap- 
proach, it is extremely unlikely that a routine VLBI observation would involve rays passing close 
enough to a gravitating body for this correction to be of importance. Exceptions are experiments 
specifically designed to measure planetary gravitational bending (TYeuhaft and Lowe, 1991). In order 
to guard against such an unlikely situation in routine work, and to provide analysis capability for spe- 
cial experiments, the MODEST code always performs the transit-time correction for all planets. To 
obtain the positions of the gravitational objects, we employ an iterative procedure, using the positions 
and velocities of the objects at signal reception time. If R(t r ) I s the position of the gravitational object 
at signal reception time, £ r , then that object’s position, R(£ a )> at the time, t ai of closest approach of 
the ray path to the object was: _ 

H{t a )=K{t r )-\[t r -ta} (2.30) 


t r t a — 



(2.31) 


We do this correction iteratively, using the velocity, V(£ r ), as an approximation of the mean velocity, 
V. Because v/c « 10" 4 , an iterative solution: 


K n {t a ) = R(t r ) - 


V(tr) 


iRn— 1(*«)| 


(2.32) 


rapidly converges to the required accuracy. 

Gravitational potential effects and curved wave front effects are calculated independently of each 
other since the gravitational effects are a small perturbation (« 8.5 microradians or < 1. 75) for 
Sun-grazing rays. 
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2.2 TIME INFORMATION 

Before continuing the description of the geometric model, a few words must be said about tune- 
tag information and the time units which will appear as arguments below. A general reference for 
time definitions is the Explanatory Supplement, 1961. The epoch timing information m the data is 
taken from the UTC (Universal Coordinated Time) time tags in the data stream at station #1. ihis 
time is converted to Terrestrial Dynamic Time (TDT) and is also used as an argument to obtain an 
a priori estimate of Earth orientation. The conversion consists of the following components: 

TDT = ( TDT - TAI) + (TAI - UTCjers) + ( UTCiers ~ UTC , ) 

+ {UTCo - UTC!) + UTCi ( 2 - 33 ) 


where in seconds: 


TDT -TAI = 32.184 


(2-34) 


and where TAI (Temps Atomique International) is atomic time. The International Earth Rotation 
Service (IERS), its predecessor, Bureau International de l’Heure (BIH), and Bureau International 
des Poids et Mesures (BIPM) are the coordinating bodies responsible for upkeep and publication of 
standard time and Earth rotation quantities. TAI -UTC IE rs = published integer second offset after 
O' 1 , January 1, 1972 (leap seconds), and 

TAI - UTCjers = 9.8922417 + 3.0 x 10 -8 x ( UTCiers ~ UTC, iers ) (2-35) 


between 0 fc , January 1, 1968, and 0\ January 1, 1972. UTCjers ~ UTC 0 iers = number of UTC 
seconds relative to January 1, 1972. This is a negative number prior to that date. The software will 
not allow this quantity to be obtained prior to 1968. UTCjers ~UTC 0 = the offset m UTC seconds 
between IERS UTC and the UTC clock at some secondary standard (usually NBS in Boulder for DSN 
observations). This can be obtained from BIPM Circular T (typical reference is Bureau International 
des Poids et Mesures, 1990). In practice as of January, 1972, all that we do is use a linear interpolation 
between (UTCifes - UTCjjbs ) data points as published in IERS Bulletin A. The approximation 
usually is made S the clock at station #1 is very close to the NBS clock, UTC, - UTCj < 
5-10 p.s. Since this time is used as epoch time in the observations, the major consequence resulting 
from an error in this assumption is to make an error in the estimation of UT1-UTC of one second per 
second of error in {UTC, - UTC } ). An error in epoch time causes an error of « = 7.3 X 10 

cm per km baseline per ps of clock error, where we is the rotation rate of the Earth (section 7). Even 
for the extreme case of a 10,000 km baseline and At = 10 /is, this amounts to only 0.007 cm. 

A priori UT1-UTC and pole positions are normally obtained by interpolation of the IERS 
Bulletin A smoothed values. However, any other source of UT1-UTC and pole position could be 
used provided it is a function of UTC, and is expressed in a left-handed coordinate system (see section 
2.6.1). Part of the documentation for any particular set of results should clearly state what were the 
values of UT1-UTC and pole position used in the data reduction process. 

For the Earth model based on the new IAU conventions, the following definitions are employed 

throughout (Kaplan, 1981): 

1. Julian date at epoch J2000 = 2451545.0. 

2. All time arguments denoted by T below are measured in Julian centuries of 36525 days of the 
appropriate time relative to the epoch J2000, i.e., T — [JD — 2451545.0)/36525. 

3. For the time arguments used to obtain precession, nutation, or to reference the ephemens, 
Barycentric Dynamic Time (TDB, Temps Dynamique Barycentrique) is used. This is related 
to Terrestrial Dynamic Time {TDT, Temps Dynamique Terrestre) by the following. 


TDB = TDT + 0. ’001658 sin(g + 0.0167 sin(g)) 


(2.36) 


where 


(357.°528 + 35999.°050 TDT) x 2x 
360° 


(2.37) 
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2.3 STATION LOCATIONS 


Coordinates of the observing stations are expressed in the Conventional International Origin 
(CIO) 1903.0 reference system, with the reference point for each antenna defined as in Sec. 2.8. The 
pre-1984 model considered the three coordinates of station i: r 9Pi , A;, Z{ (radius off spin axis, longitude, 
and height above the equator, respectively) to be time- invariant. In investigations of tectonic motion, 
however, a new set of coordinates is usually solved for in the least-squares estimation process for each 
VLBI session. Post-processing software then makes linear fits to these results to infer the time rate of 
change of the station location. Care must be taken that the correlations of coordinates estimated at 
different epochs are accounted for properly. The advantage of this approach is that the contribution of 
each session to the overall slope may be independently evaluated, since it is clearly isolated. Since this 
procedure is somewhat inconvenient in practice, an alternative is to introduce the time rates of change 
of the station coordinates as new parameters in MODEST. The model is linear, with the cylindrical 
coordinates at time t expressed as 


r »Pi — r ap , + t 0 ) (2.38) 

A i = A 9 + A<(t - t 0 ) (2.39) 

= z? +Zi(t - t 0 ) (2.40) 

Here to is a reference epoch, at which the station coordinates are A9, #9). If modeling is done 

in Cartesian coordinates, the analogous expressions are 

Xi = x? + Xi(t- to) (2.41) 

!/i = y° + in(t-t 0 ) (2.42) 

= 2° + *i(t - <o) (2.43) 

with (x°, t/9 , z?) being the station coordinates at the reference epoch. 

2.3.1 Models of Tectonic Plate Motion 

As an alternative to estimating linear time dependence of the station coordinates, two standard 
models of tectonic plate rotation are optionally available in MODEST. The first is described in an 
addition to the MERIT standards document (Melbourne et ah, 1985), and was denoted AMO-2 in the 
original paper (Minster and Jordan, 1978). Time dependence of the Cartesian station coordinates is 
expressed as 


Xi = I t ° + - f 0 ) 

(2.44) 

y i = y . 9 + (“i *. 9 - uiz?){t - M 

(2.45) 

z, = Z° + (u>iy° - - i 0 ) 

(2.46) 


where u J x y z are velocities of the plate j on which station i resides. Table I gives a list of the rotation 
rates for the 11 plates in the AMO-2 model. 
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Table I 

Plate Rotation Velocities: Minster- Jordan AMO-2 Model^ 


Plate 


Wy 

w z 

AFRC 

0.988 

-3.360 

4.192 

ANTA 

-0.923 

-1.657 

3.765 

ARAB 

4.867 

-2.922 

6.520 

CARB 

-0.486 

-0.988 

1.881 

COCO 

-11.122 

-23.238 

12.663 

EURA 

-0.536 

-2.769 

3.422 

INDI 

8.443 

4.365 

7.528 

NAZC 

-1.586 

-9.299 

11.006 

NOAM 

0.576 

-3.984 

-0.249 

PCFC 

-2.143 

5.439 

-11.438 

SOAM 

-0.978 

-1.863 

-1.508 


t units are nrad/year 


Note that the velocities are expressed in nanoradians per year rather than the microdegrees per year 
used in the original paper. 

More recent models, denoted NUVEL-1 and NNR-NUVEL1, are due to DeMets et al. (1990) 
and Argus and Gordon (1991), respectively. In NUVEL-1, the Pacific plate is stationary, while 
NNR-NUVEL1 is based on the imposition of a no-net-rotation (NNR) condition. With some notable 
exceptions, the NUVEL models give rates that are very close to those of the AMO-2 model. The 
AMO-2 INDI plate has been split into AUST and INDI, and there are two additional plates: JDEF (Juan 
de Fuca) and PHIL (Philippine). The NUVEL-1 rotation rates are given in Tables II and III. 


Table II 

Plate Rotation Velocities: NUVEL-1 Model t 


Plate 

w * 


w* 

AFRC 

2.511 

-8.303 

14.529 

ANTA 

0.721 

-6.841 

14.302 

ARAB 

8.570 

-5.607 

17.496 

AUST 

9.777 

0.297 

16.997 

CARB 

1.393 

-8.602 

12.080 

COCO 

-9.323 

-27.657 

21.853 

EURA 

0.553 

-7.567 

13.724 

INDI 

8.555 

-5.020 

17.528 

JDEF 

6.81 

3.32 

5.31 

NAZC 

-0.023 

-14.032 

20.476 

NOAM 

1.849 

-8.826 

10.267 

PCFC 

0.000 

0.000 

0.000 

PHIL 

11.9 

12.8 

0.000 

SOAM 

0.494 

-6.646 

9.517 


| units are nrad/year 
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Table III 

Plate Rotation Velocities: NNR-NUVELl Modelt 



f units are nrad/year 


At present, there is no facility in MODEST to compute partial derivatives with respect to the 
plate velocities, or to solve for these quantities. 
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2.4 TIDAL EFFECTS 


As an initial step in calculating the geometric delay, we need to consider the effects of crustal 
motions on station locations. Among these deformations are solid Earth tides, tectonic motions, and 
alterations of the Earth’s surface due to local geological, hydrological, and atmospheric processes. One 
possibility is to not model crustal movement other than that due to solid Earth tides, and allow the 
other effects to manifest themselves as temporal changes of the Earth-fixed baseline. Such a strategy 
corrupts the estimation of global orientation parameters from a finite set of baselines, and is a known 
weakness (wl-10 cm/year) of the simplified form of the current model. 

In the standard terrestrial coordinate system, tidal effects modify the station location ro by an 
amount 

A = A soi + -h (2.47) 

where the four terms are due to solid Earth tides, pole tide, ocean loading, and atmosphere loading, 
respectively. Other Earth-fixed effects would be incorporated by augmenting the definition of A. 
All four tidal effects are most easily calculated in some variant of the VEN (vertical, East, North) 
local geocentric coordinate system. To transform them to the Earth-fixed coordinate frame, the 
transformation VW f given in the next section, is applied. 

2.4.1 Solid Earth Tides 

Calculating the alteration of the positions of the stations caused by solid Earth tides is rather 
complicated due to the solid tides’ coupling with the ocean tides, and the effects of local geology. 
We have chosen to gloss over these complications initially, and to incorporate the simple multipole 
response model described by Williams (1970), who used Melchior (1966) as a reference. Let R p be the 
position of a perturbing source in the terrestrial reference system, and r 0 the station position in the 
same coordinate system. To allow for a phase shift (ip) of the tidal effects, the phase-shifted station 
vector t 9 is calculated from tq by applying a right-handed rotation, L, through an angle ip about the 
Z axis of date, r 9 = Lro. This lag matrix, L } is: 

( cos ip sin ip 0 \ 

— sin ip cos ip 0 J (2.48) 

ooi ; 

By a positive value of ip we mean that the peak response on an Earth meridian occurs at a time 
St = ip /we after that meridian containing r 0 crosses the tide-producing object, where we is the 
angular rotation rate of the Earth. In the vertical component, the peak response occurs when the 
meridian containing t 9 also includes R p , 

The tidal potential at r, due to the perturbing source at R p is expressed as 

= U 2 + U 3 (2.49) 

where only the quadrupole and octupole terms have been retained. Here, G is the gravitational 
constant, m p is the mass of the perturbing source, Pi are the Legendre polynomials, and 9 is the angle 
between t 9 and R p . 

In a local geocentric VEN coordinate system (axes vertical, eastward, and northward) on a 
spherical Earth, the tidal displacement vector 6 is 

* = £>f (2-50) 
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where the (t = 2, 3) are the quadrupole and octupole displacements. The components of S are 
obtained from the tidal potential as 


gf = kiUi/g 








(2.51) 

(2.52) 

(2.53) 


where h, (i = 2,3) are the vertical (quadrupole and octupole) Love numbers, I, (t = 2,3) the corre- 
sponding horizontal Love numbers, and X 8 and <f> 9 are the station longitude and latitude, and g the 
acceleration due to gravity, 

g = Gm E /r * (2-54) 

Using the relation between terrestrial and celestial coordinates, 


cos 0 — sin <f> 8 sin S p + cos <f> a cos 6 p cos (A, + cxq — a p ) (2.55) 

with a p , S p the right ascension and declination of the perturbing body, and c*g the RA of Greenwich, 
some algebra produces the following expressions for the quadrupole and octupole components of S in 
terms of the coordinates of the station (x Sl y 3i z s ) and the tide-producing bodies (X p , Y P) Z p ): 


-(2) _ Y' 3 /b» r ? r ( r » ■ R p f r > R P 1 

31 "^“rTL 2 6 J 


s4 2) = E ^( r * *p)1*>Yp - 


4 2) - E ( r * R p) W x * + y? Z P - — J X > X P + y> y p ) 


\Z^f+yf 




5(r, ■ Rp) 2 — 3r^fip 


J3) _ 3 ^P r « 

^2 2R 7 

p Zn P L 


5(r» * Rp) 2 — 


(*,y p - y .x p )AAfT^F 


4 3) - E [ 5 ( r * • R p) 2 - [\^T+^ z p - [x.x p + y.Y p ) 

p zn p L j L V 1 ? + y; 


(2.56) 

(2.57) 

(2.58) 

(2.59) 

(2.60) 
(2.61) 


where fi p is the ratio of the mass of the disturbing object, p, to the mass of the Earth, and 

R p ~ [Xp, Y p , Z p ] (2.62) 

is the vector from the center of the Earth to that body. The summations are over tide-producing 
bodies, of which we include only the Sun and the Moon. If the tidal effect at time ti is desired, and 
the light travel time is St, then the position of the tide- producing mass at time 

ti - St = ti - |R p (tx - St)\/c (2.63) 

should be used (a nuance we have not yet incorporated). While the quadrupole displacements are 
of the order of 50 cm, the mass and distance ratios of the Earth, Moon, and Sun limit the octupole 
terms to a few mm. The octupole terms are optionally included in the MODEST code, but partials 
with respect to the Love numbers are available only for the quadrupole terms. 
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To convert the locally referenced strain, 6, which is expressed in the VEN system, to the Earth- 
fixed frame, two rotations must be performed. The first, W , rotates by an angle, <f>» (station geodetic 
latitude), about the y axis to an equatorial system. The second, V , rotates about the resultant z axis 
by angle, -A, (station longitude), to bring the displacements into the standard geocentric coordinate 
system. The result is 

A <0 , = VW6 (2.64) 


where 


and 



0 — sin <f> 9 \ 

1 0 

0 cos <f> 9 J 



(2.65) 


( 2 . 66 ) 


Actually, the product of these two matrices is coded: 


VW 


( cos A, cos <f> 9 
sin A* cos <f> 8 
sin<^ 


— sin A, 
cos A, 
0 


COS ( 


(2.67) 


MODEST code uses geodetic latitudes 


( 2 . 68 ) 


where / is the geoid flattening factor. The difference between geodetic and geocentric latitude can 
affect this model on the order of (tidal effect) /(flattening factor) » 0.1 cm. 


2.4.2 Pole Tide 


One of the secondary tidal effects is the displacement of a station by the elastic response of the 
Earth’s crust to shifts in the spin axis orientation. The spin axis is known to describe a circle of 
« 20-m diameter at the north pole. Depending on where the spin axis pierces the crust at the instant 
of a VLBI measurement, the “pole tide” displacement will vary from time to time. This effect must 
be included if centimeter accuracy is desired. 

Yoder (1984) derived an expression for the displacement of a point at geocentric latitude <f> t 
longitude A due to the pole tide: 


a/ 2 R 

6 = — — [sin <f> cos <f>(x cos A 4* y sin A) hr 

+ cos 2 <f>(x cos A *f y sin A) If 

+ sin <£(— x sin A -h y cos A) l A] (2.69) 


Here we is the rotation rate of the Earth, R the radius of the (spherical) Earth, g the acceleration 
due to gravity at the Earth’s surface, and h and l the customary Love numbers. Displacements of the 
spin axis from the 1903.0 CIO pole position along the x and y axes are given by x and y. Eq. (2.69) 
shows how these map into station displacements along the unit vectors in the radial (r), latitude (^), 
and longitude (A) directions. With the standard values we ~ 7.292 X 10 -5 rad/sec, R = 6378 km, 
and g = 980.665 cm/sec 2 , the factor w\Rjg = 3.459 X 10" 3 . Since the maximum values of x and y 
are of the order of 10 meters, and h » 0.6, / w 0.08, the maximum displacement due to the pole tide 
is 1 to 2 cm, depending on the location of the station (<{>, A). 

The locally referenced displacement 6 is transformed via the suitably modified transformation 
(2.67) to give the displacement A po i in the standard geocentric coordinate system. The pole tide 
effect has been coded as an optional part of the MODEST model. It is only applied if specifically 
requested, t.e., the default model contains no pole tide contributions to the station locations. 
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2.4.3 Ocean Loading 

This section is concerned with another of the secondary tidal effects, i.e., the elastic response 
of the Earth’s crust to ocean tides, which move the observing stations to the extent of a few cm. 
Such effects are commonly labeled “ocean loading.” A model of ocean loading is incorporated in 
the MODEST code. It is general enough to accommodate a variety of externally derived constants 
describing the tide phases and amplitudes. Because the station motions caused by response to ocean 
tides appear to be limited to approximately 3 cm for sites further than «100 km from the coast, no 
estimation capability was deemed necessary at present. This decision is supported by the fact that for 
locations near the coast, where the effects may be more sizeable, and which would thus be expected 
to produce data useful in parameter estimation, the elastic response modeling is as yet inadequate 
(Agnew, 1982). As suggested in section 8 of the initial version of this report (Fanselow, 1983), local 
Earth motion can be partially accounted for by varying the Love numbers for each station. The 
present model entails deriving an expression for the locally referenced displacement S due to ocean 
loading. In the vertical, N-S, E-W local coordinate system (the computer code accepts inputs related 
to unit vectors in the vertical, North, and West directions) at time t, 

N 

Sj = y ^ g? cos(fa>,-t + Vj- S-) (2.70) 

« = 1 

The quantities a;* (frequency of tidal constituent i) and V{ (astronomical argument of constituent t) 
depend only on the ephemeris information (positions of the Sun and Moon). The algorithm of Goad 
(IERS, 1989) is used to calculate these two quantities. On the other hand the amplitude £? and 
Greenwich phase lag 6* of each tidal component j are determined by the particular model assumed 
for the deformation of the Earth. The local displacement vector is transformed via Eqs. (2.67) and 
(2.64) to the displacement A ocn in the standard geocentric frame. 

Input to MODEST provides for specification of up to 11 frequencies and astronomical arguments 
oji and V t , followed by tables of the local distortions and their phases, £? and £?, calculated from the 
ocean tidal loading model of choice. The eleven components are denoted, in standard notation: Af 2 , 
S 2) N 2} and K 2 (all with approximately 12-hour periods), K X} 0 X) P x , Q x (24 hr), Mj (14 day), Af m 
(monthly), and S 8a (semiannual). 

Presently four choices of ocean loading models are available for use with MODEST. They differ in 
the displacements calculated and components considered, as well as in the numerical values that they 
yield for the £?s and Sf s. Schemeck’s results (1983, 1990, 1991) are the most complete in the sense 
of considering both vertical and horizontal displacements and all eleven tidal components. Goad’s 
model (1983) has been adopted in the MERIT and IERS standards (1989), but only considers vertical 
displacements. Pagiatakis’ (1982, 1990) model, based on Pagiatakis, Langley, and Vanicek (1982), 
considers only six tidal components. Agnew (1982) only considers five components, but pays special 
attention to points near coastlines. Table IV summarizes the features of the four models, with V and 
H indicating vertical and horizontal components, respectively. 

Due to their bulk, none of the tables of tidal amplitudes is reproduced here, but are available 
on request in computer-readable form. The default tidal model in MODEST remains the Williams 
quadrupole solid Earth tide model with no ocean loading. 


Table IV. Ocean Loading Models 


Model 

Displacements 

Tidal components 

Scherneck 

V, H 

M 2 S 2 N 2 K 2 KyO l P l Q l M J M rn S ia 

Goad (MERIT, IERS) 

V 

M 2 S 2 N 2 K 2 K\0\PiQ\MfM rn S 9a 

Pagiatakis 

V, H 

M 2 S 2 N 2 K x O x Pi 

Agnew 

V, H 

M 2 S 2 N 2 K 1 O l 
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2.4.4 Atmosphere Loading 


By analogy with the consequences of ocean tides that were considered in the previous section, a 
time-varying atmospheric pressure distribution can induce crustal deformation. A paper by Rabbel 
and Schuh (1986) estimates the effects of atmospheric loading on VLBI baseline determinations, and 
concludes that they may amount to many millimeters of seasonal variation. In contrast to ocean 
tidal effects, analysis of the situation in the atmospheric case does not benefit from the presence of 
a well-understood periodic driving force. Otherwise, estimation of atmospheric loading via Green’s 
function techniques is analogous to methods used to calculate ocean loading effects. Rabbel and 
Schuh recommend a simplified form of the dependence of the vertical crust displacement on pressure 
distribution. It involves only the instantaneous pressure at the site in question, and an average pressure 
over a circular region C of radius R = 2000 km surrounding the site. The expression for the vertical 


displacement (mm) is: 


Ar = — O.35po — 0.55p 


(2.71) 


where po is the local pressure anomaly (relative to the standard pressure of 1013 mbar), and p 
the pressure anomaly within the 2000-km circular region mentioned above (both quantities are in 
mbar). Note that the reference point for this displacement is the site location at standard (sea level) 
pressure. The locally referenced Ar is transformed to the standard geocentric coordinate system via 
the transformation (2.67). 

It was decided to incorporate this rudimentary model into MODEST as an optional part of the 
model, with an additional mechanism for characterizing p. The two-dimensional surface pressure 
distribution (relative to 1013 mbar) surrounding a site is described by 


p(x , y) = Aq + A\x + A 2 y + Asx 2 + A^xy + A$y 2 


(2.72) 


where x and y are the local East and North distances of the point in question from the VLBI site. 
The pressure anomaly p may then be evaluated by the simple integration 


= JJ dx dy p(x, y) / J J c 


p — Aq 4- (A3 + A&)R 2 /4 (2-74) 

It remains the task of the data analyst to perform a quadratic fit to the available weather data 
to determine the coefficients A0-5. Future advances in understanding the atmosphere-crust elastic 
interaction can probably be accommodated by adjusting the coefficients in Eq. (2.71). 

After each of the locally referenced tidal displacements has been transformed to standard terres- 
trial coordinates, the station location is 

r t = To + 4* A po | + Aocn H" A atm (2.75) 
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2.5 SOURCE STRUCTURE EFFECTS 


Numerous as trophy sic al studies during the past decade have shown that compact extragalactic 
radio sources have structures on a milliarcsecond scale (e.g., Kellermann and Pauliny-Toth, 1981), 
Such studies are important for developing models of the origin of radio emission of these objects. 
Many radio source structures are found to be quite variable with frequency and time (Zensus and 
Pearson, 1987). If extragalactic sources are to serve as reference points in a stable reference frame, it 
is important to correct for the effects of their structures in astrometric VLBI observations. 

Recently, MODEST modeling was extended to allow optional corrections for the effects of source 
internal structures, based on work by Thomas (1980), Ulvestad (1988), and Chariot (1989). A non- 
point like distribution of the intensity of a source yields time dependent corrections to the group delay 
and delay rate observables, A r 9 and A r 91 that may be written in terms of the intensity distribution 
J(s,a;,t) as 

Ar 9 = d<f> 9 /dw , A r 9 =d<f> 9 /dt (2.76) 

with 

<f>s = arctan(— 2^/Zc) (2.77) 

and 

Z { ‘}= J j dn 7(8, o>, t) { ™ } ( 2 ttB s/A) (2.78) 

Here <f> 9 is the correction to the phase of the incoming signal, s is a vector from the adopted reference 
point to a point within the source intensity distribution in the plane of the sky, w and A are the 
observing frequency and wavelength, B the baseline vector, and the integration is over solid angles D. 
Source intensity distribution maps are most conveniently parametrized in terms of one of two models: 
superpositions of delta functions or Gaussians. At a given frequency, the corresponding intensity 
distributions are written as 

-^(s) = — v — VJb) (2.79) 

k 

or 


O r 

7 ( 8 ) = 2na k b k exp l K* ~ Tfc ) cos + (y ~ y/t) sin 6 k ] 2 /2a k 2 

-(-(* — Sfc) sin 6 k + (y - y k ) cos 0 k \ 2 / 26 fc 2 j ( 2 . 80 ) 

where S k is the flux of component k } and (with components y * in the plane of the sky) is its 

position relative to the reference point. For Gaussian distributions, 6 k is the angle between the major 
axis of component k and the u axis (to be defined below), and (a*, 6*) are the full widths at half 
maximum of the (major, minor) axes of component k normalized by 2^/2~Iog~2. The quantities -£{'} 
entering the structure phase <j> 9 [Eq. (2.77)] are 

: ) = Y, s « { “ } ( 2 * B • s “/ A ) (2- 81 ) 

k 


for delta functions, and 




-5> 


k ex P[ 


-2?r' 


’■(4u k 


2 + W 


)1{™}(2, B ../A) 


(2.82) 


for Gaussians. Here 

Uk = u cos Ok -f- usin^fc (2.83) 

V k = -usinO k + v cos 6 k (2. 84) 

with u, v being the projections of the baseline vector B on the plane of the sky in the E-W, N-S 
directions, respectively. 
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MODEST accepts maps specified in terms of an arbitrary number of Gaussian or delta function 
components. At most, six parameters must be specified for each component: its polar coordinates 
and flux, and, for a Gaussian, its major and minor axes and the position angle of the major axis. The 
structural correction for phase is computed via Eqs. (2.77), (2.81), and (2.82). For the BWS delay 
observable, the structure correction is the slope of a straight line fitted to the individual structure 
phases calculated for each frequency channel used during the observation. For example, for Mark III 
data there are typically 8 channels spanning «8.2 to 8.6 GHz at X band, and 6 channels spanning 
«2.2 to 2.3 GHz at S band. Delay rate structure corrections are calculated by differencing the 
structure phases at i2 seconds (see Section 6). In the case of dual-band (S-X) experiments, a linear 
combination of the structure corrections calculated independently for each band is applied to the 
dual-band observables. 

The practical question to be resolved in the future is whether such structural corrections yield 
significant and detectable corrections to the observables at the present levels of experimental and 
modeling uncertainty. Maps are available for only a few of the hundreds of sources currently observed 
by VLBI. Some of the extended sources show time variability on a scale of months; since the corrections 

A t s and A t 9 are quite sensitive to fine details of the structure, in such cases new maps may be 

required on short time scales. Depending on the relative orientation of the source and baseline, the 
delay correction can be as large as ^1 ns, which is equivalent to tens of cm. An optimistic note is the 
recent observation of Chariot (1990) that data from a multiple baseline geodynamics experiment are 
adequate to map source structures with high angular resolution. 

Empirical evaluation of the effects of unknown source structure on VLBI measurements could be 
made via the time rates of change of the source right ascension a and declination 6, A linear model 
of the motion of source coordinates 

a = olq H- a(t — to) (2.85) 

6 = 6 0 + 6 {* - *o) (2-86) 

is implemented in MODEST. Non-zero estimates of the rate parameters a and 6 could arise either 
from genuine proper motion or from motion of the effective source centroid sampled by VLBI mea- 
surements. Proper interpretation of such results is problematic, but non-zero rates can be used as a 
crude diagnostic for the presence of structure effects. 
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2.6 TRANSFORMATION FROM TERRESTRIAL TO CELESTIAL 
COORDINATE SYSTEMS 

The Earth is approximately an oblate spheroid, spinning in the presence of two massive moving 
objects (the Sun and the Moon) which are positioned such that their time-varying gravitational effects 
not only produce tides on the Earth, but also subject it to torques. In addition, the Earth is covered 
by a complicated fluid layer, and also is not perfectly rigid internally. As a result, the orientation 
of the Earth is a very complicated function of time, which to first order can be represented as the 
composite of a time-varying rotation rate, a wobble, a nutation, and a precession. The exchange of 
angular momentum between the solid Earth and the fluids on its surface is not readily predictable, 
and thus must be continually determined experimentally. Nutation and precession are well modeled 
theoretically. However, at the accuracy with which VLBI can determine baseline vectors, even these 
models are not completely adequate. 

Currently, the rotational transformation, Q , of coordinate frames from the terrestrial frame to 
the celestial geocentric frame is composed of 6 separate rotations (actually 12, since the nutation, 
precession, and “perturbation” transformations, N> P, and fi, consist of 3 transformations each) 
applied to a vector in the terrestrial system: 

Q - OPNUXY (2.87) 

In order of appearance in (2.87), the transformations are: the perturbation rotation, precession, 
nutation, UT1, and the x and y components of polar motion. All are discussed in detail in the 
following four sections. With this definition of Q> if r« is a station location expressed in the terrestrial 
system, e.g., the result of (2.75), that location, r c , expressed in the celestial system is 

r c = Qr t (2.88) 

This particular formulation follows the historical path of astrometry, and is couched in that 
language. While esthetically unsatisfactory with modern measurement techniques, such a formulation 
is currently practical for intercomparison of techniques and for effecting a smooth inclusion of the 
interferometer data into the long historical record of astrometric data. Much more pleasing esthetically 
would be the separation of Q into two rotation matrices: 

Q - Q 1 Q 2 (2.89) 

where Q 2 are those rotations to which the Earth would be subjected if all external torques were 
removed (approximately UXY above), and where Q 1 are those rotations arising from external torques 
(approximately QPN above). Even then, the tidal response of the Earth prevents such a separation 
from being perfectly realized. Eventually, the entire problem of obtaining the matrix Q , and the tidal 
effects on station locations should be done numerically. Note that the six rotations operating on a 
vector yield its components in a new coordinate system, and, since we rotate the Earth rather than the 
celestial sphere, the matrices H, P, and N will be the transposes of those used to rotate the celestial 
system of J2000 to a celestial system of date. 

2.6.1 UT1 AND POLAR MOTION 

The first transformation, Y, is a right-handed rotation about the x axis of the terrestrial frame 
by an angle 02- Currently, the terrestrial frame is the 1903.0 CIO frame, except that the positive y 
axis is at 90 degrees east (Moscow). The x axis is coincident with the 1903.0 meridian of Greenwich, 
and the z axis is the 1903.0 standard pole. 

0 "1 

sin 0 2 I (2.90) 

cos 02 J 


1 0 

Y = I 0 cos 02 
l 0 — sin 02 
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where ©2 is the y pole position published by IERS. 

The next rotation in sequence is the right-handed rotation through an angle ©1 about the y axis 

obtained after the previous rotation has been applied: 


( cos©i 0 

0 i 

sin©i 0 


— sin ©i 

0 

COS ©1 


(2.91) 


In this rotation, ©i is the IERS x pole position. Note that we have incorporated in the matrix 
definitions the transformation from the left-handed system nsed by IERS to the right-handed system 
we use. Note also that instead of IERS data used as a pole definition, we could instead use any other 
source of polar motion data provided it was represented in a left-handed system. The only effect 
would be a change in the definition of the terrestrial reference system. 

The application of U XY” to a vector in the terrestrial system of coordinates expresses that vector 
as it would be observed in a coordinate frame whose z axis was along the Earth’s ephemeris pole. 
The third rotation, U, is about the resultant z axis obtained by applying *XY*. It is a rotation 
through the angle, -H, where H is the hour angle of the true equinox of date (».e., the dihedral angle 
measured westward between the xz plane defined above and the meridian plane containing the true 
equinox of date). The equinox of date is the point defined on the celestial equator by the intersection 
of the mean ecliptic with that equator. It is that intersection where the mean ecliptic rises from below 
the equator to above it (ascending node). 


( COB H 
sin H 
0 


— sin H 
cos H 
0 



This angle H is composed of two parts: 


H — + as 


(2.92) 


(2.93) 


where /i 7 is the hour angle of the mean equinox of date, and cie (equation of equinoxes) is the difference 
in hour angle of the true equinox of date and the mean equinox of date, a difference which is due 
to the nutation of the Earth. This set of definitions is cumbersome and couples the nutation and 
precession effects into Earth rotation measurements. However, in order to provide a direct estimate 
of conventional UT 1 it is convenient to endure this historical approach, at least for the near future.. 

UTl (universal time) is defined to be such that the hour angle of the mean equinox of date is 
given by the following expression (Aoki et al. t 1982, and Kaplan, 1981): 

h y = UTl + 6 h 41 m 50*. 54841 + 8640184* .812866 T„ 

+ 0* .093 104 T* - 6* .2 X 10 -6 (2.94) 


where the dimensionless quantity 

_ (Julian UTl date) - 2451545.0 
Tu ~ 36525 


(2.95) 


The actual equivalent expression which is coded is: 

—2k{UT1 Julian day fraction) + 67310®. 54841 

+ 8640184 s . 812866 T u + 0 S .093104 - 6\2 X 10” 6 Tj (2.96) 

This expression produces a time, UTl, which tracks the Greenwich hour angle of the real Sun to 
within 16 m . However, it really is sidereal time, modified to fit our intuitive desire to have the Sun 
directly overhead at noon on the Greenwich meridian. Historically, differences of UTl from a uniform 
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measure of time, such as atomic time, have been used in specifying the orientation of the Earth. Note 
that this definition has buried in it the precession constant since it refers to the mean equinox of date. 

By the very definition of “mean of date” and “true of date”, nutation causes a difference in the 
hour angles of the mean equinox of date and the true equinox of date. This difference, called the 
“equation of equinoxes”, is denoted by as and is obtained accordingly: 


where the vector 



(2.97) 


(2.98) 


is the unit vector, in true equatorial coordinates of date, toward the mean equinox of date. In mean 
equatorial coordinates of date, this same unit vector is just (1,0, 0) r . The matrix N~ l is just the 
inverse (or equally, the transpose) of the transformation matrix N which will be defined below [Eq. 
(2.105)] to effect the transformation from true equatorial coordinates of date to mean equatorial 
coordinates of date. 


2. 6. 1.1 Short Period UTl Variations 


Depending on the smoothing used to produce the a priori UTl — UTC series, the short-period 
(t < 35 days) fluctuations in UTl due to changes in the latitude and size of the mean tidal bulge may 
or may not be smoothed out. Since we want as accurate an a priori as possible, it may be necessary 
to add this effect to the UTl a priori obtained from the series UTl 9tnoothed . If this option is selected, 
then the desired a priori UTl is given by 


UTl a priori — UTl smoothed + A£7TT (2.99) 

UTl moo thtd represents an appropriately smoothed a priori measurement of the orientation of the 
Earth (i.c., typically IERS Bulletin A smoothed or, even better, *7Tlf*), for which the short period 
(t < 35 days) tidal effects have either been averaged to zero, or, as in the case of UTlR, removed 
before smoothing. This AUTl can be represented as 


N 

AUTl = ^2 

i=l 


Ai sin 


f 5 

5Z k U a i 


L y=i 


( 2 . 100 ) 


where N is chosen to include all terms with a period less than 35 days. There are no other con- 
tributions until a period of 90 day3 is reached. However, these long-period terms are included by 
the measurements of the current Earth-orientation measurement services. The values for fc,-, and A , , 
along with the period involved, are given in Table V. The a t for * = 1, 5 are just the angles defined 
below (Section 2.6.2) in the nutation series as l, l', F, D, and 0, respectively. In Table V, the sign 
of the 14.73 day term has been changed [Yoder (1982)] to correct a sign error in Yoder et al. (1981). 
The BIH Annual Report for 1982 [BIH (1983)] is the first reference to give the correct table. 

It might be appropriate at this point to describe the interpolation method used in MODEST to 
obtain a priori polar motion and UTl values. These are normally available as tables at 5-day intervals, 
from either IERS (IERS, 1991) or the IRIS project (IAG, 1986). Linear interpolation is performed 
for all three quantities. If the short-period tidal terms AUTl are present in the tabular values, they 
are subtracted before interpolation, and added back to the final value. ^Yith the present accuracy of 
determinations of pole position and UTl (1 mas and 0.05 ms respectively), linear interpolation over a 
5-day interval may be inadequate, possibly giving rise to 0.1 to 0.2 ms errors in UTl. Quadratic spline 
interpolation is being considered as an alternative. Even with the present code, however, the highest 
possible accuracy may be achieved by performing the interpolation externally to MODEST, and 
supplying it with tables of values more closely spaced in time for the final internal linear interpolation. 
The Kalman-filtered UTPM values of Eubanks et al. (1984) are ideally suited for this purpose. 
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Table V 

Periodic Tidally Induced Variations in UT1 
with Periods Less than 35 Days 


Index 

i 

Period 

(days) 

Argument coefficient 

&i'2 ^*3 4 5 

Ai 

(O'.OOOl) 

1 

5.64 

1 

0 

2 

2 

2 

-0.02 

2 

6.85 

2 

0 

2 

0 

1 

-0.04 

3 

6.86 

2 

0 

2 

0 

2 

-0.10 

4 

7.09 

0 

0 

2 

2 

1 

-0.05 

5 

7.10 

0 

0 

2 

2 

2 

-0.12 

6 

9.11 

1 

0 

2 

0 

0 

-0.04 

7 

9.12 

1 

0 

2 

0 

1 

-0.41 

8 

9.13 

1 

0 

2 

0 

2 

-0.99 

9 

9.18 

3 

0 

0 

0 

0 

-0,02 

10 

9.54 

-1 

0 

2 

2 

1 

-0.08 

11 

9.56 

-1 

0 

2 

2 

2 

-0.20 

12 

9.61 

1 

0 

0 

2 

0 

-0.08 

13 

12.81 

2 

0 

2 

-2 

2 

0.02 

14 

13.17 

0 

1 

2 

0 

2 

0.03 

15 

13.61 

0 

0 

2 

0 

0 

-0.30 

16 

13.63 

0 

0 

2 

0 

1 

-3.21 

17 

13.66 

0 

0 

2 

0 

2 

-7.76 

18 

13.75 

2 

0 

0 

0 

-1 

0.02 

19 

13.78 

2 

0 

0 

0 

0 

-0.34 

20 

13.81 

2 

0 

0 

0 

1 

0.02 

21 

14.19 

0 

-1 

2 

0 

2 

-0.02 

22 

14.73 

0 

0 

0 

2 

-1 

0.05 

23 

14.77 

0 

0 

0 

2 

0 

-0.73 

24 

14.80 

0 

0 

0 

2 

1 

-0.05 

25 

15.39 

0 

-1 

0 

2 

0 

-0.05 

26 

23.86 

1 

0 

2 

-2 

1 

0.05 

27 

23.94 

1 

0 

2 

-2 

2 

0.10 

28 

25.62 

1 

1 

0 

0 

0 

0.04 

29 

26.88 

-1 

0 

2 

0 

0 

0.05 

30 

26.98 

-1 

0 

2 

0 

1 

0.18 

31 

27.09 

-1 

0 

2 

0 

2 

0.44 

32 

27.44 

1 

0 

0 

0 

-1 

0.53 

33 

27.56 

1 

0 

0 

0 

0 

-8.26 

34 

27.67 

1 

0 

0 

0 

1 

0.54 

35 

29.53 

0 

0 

0 

1 

0 

0.05 

36 

29.80 

1 

-1 

0 

0 

0 

-0.06 

37 

31.66 

-1 

0 

0 

2 

-1 

0.12 

38 

31.81 

-1 

0 

0 

2 

0 

-1.82 

39 

31.96 

-1 

0 

0 

2 

1 

0.13 

40 

32.61 

1 

0 

-2 

2 

-1 

0.02 

41 

34.85 

-1 

-1 

0 

2 

0 

-0.09 
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It is convenient to apply U UXY ” as a group. To parts in 10 12 , XY = YX. However, with the 
same accuracy UXY / XYU. Neglecting terms of 0(© 2 ) (which produce station location errors of 
approximately 6 x 10" 4 cm): 


UXY = 


cos H 
sin H 
sin ©i 


— sin H 
cos H 

— sin ©2 


sin © i cos . 


( 2 . 101 ) 


2 . 6.2 NUTATION 

With the completion of the UTl and polar motion transformations, we are left with a station 
location vector, T da te- This is the station location relative to true equatorial celestial coordinates of 
date. The last set of transformations are nutation, N } precession, P, and the perturbation rotation, 
fl, applied in that order. These transformations give the station location, r c , in celestial equatorial 
coordinates: 

t c = nPNr daU (2.102) 

The transformation matrix N is a composite of three separate rotations (Melbourne ct a/., 1968): 

1. A(e): true equatorial coordinates of date to ecliptic coordinates of date. 

/ 1 0 M 

A(e ) = I 0 cose sine I (2.103) 

\0 — sin e cos e J 

2. C T (Sip): nutation in longitude from ecliptic coordinates of date to mean ecliptic coordinates of 
date. 


cos Sip sin Sip 0 
C T (Sip) = J —sin Sip cos Sip 0 
0 0 1 

where Sip is the nutation in ecliptic longitude. 


(2.104) 


3. A*(e) : ecliptic coordinates of date to mean equatorial coordinates. 

In ecliptic coordinates of date, the mean equinox is at an angle Sip = Ss — e — e 

is the nutation in obliquity, and c is the mean obliquity (the dihedral angle between the plane of the 
ecliptic and the mean plane of the equator). “Mean” as used in this section implies that the short- 
period ( T < 18.6 years) effects of nutation have been removed. Actually, the separation between 
nutation and precession is rather arbitrary, but historical. The composite rotation is: 

N = A T {e)C T {6i/>)A{e) (2.105) 

( cos Sip cos 5 si nSip sins sin Sip 

— cos?sin5i/> cos « cos £ cos Sip *+■ sine sine cos esin ecos Sip — sine cose 

— ainesinSip sin ecos ecos S ip — cosesine sin e sin e cos S ip + cosecose 

The 1980 IAU nutation model (Seidelmann, 1982, and Kaplan, 1981) is used to obtain the values 
for Sip and e — e. The mean obliquity is obtained from Lieske et al. (1977) or from Kaplan (1981): 


e = 23° 26' 21 "448 - 46 "8150 T - 5."9 x 10~ 4 T 2 + l."813 x 10“ 3 T 3 


T = 


(Julian TDB date) - 2451545.0 
36525 


(2.106) 

(2.107) 
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This nutation in longitude (Sip) and in obliquity ( 6s = s-s ) can be represented by a series expansion 
of the sines and cosines of linear combinations of five fundamental arguments. These are (Kaplan, 
1981, Cannon, 1981): 

1. the mean anomaly of the Moon: 

ai = l= 485866”. 733 + (1325 r + 715922".633) T 

+ 31". 310 T 2 + 0".064 T 3 (2.108) 


2. the mean anomaly of the Sun: 

a 2 = l' = 1287099". 804 + (99 r + 129258l".224) T 

- 0".577 T 2 - 0 "012 T 3 

3. the mean argument of latitude of the Moon: 

a 3 = F = 335778". 877 + (1342 r + 295263".137) T 

- 13".257 T 2 + 0".011 T 3 


(2.109) 


( 2 . 110 ) 


4. the mean elongation of the Moon from the Sun: 

a 4 = D= 1072261". 307 + (1236 r + 110560l".328) T 
- 6". 891 T 2 + 0".019 T 3 


( 2 . 111 ) 


5. the mean longitude of the ascending lunar node: 

a 6 = n = 450160".280 - (5 r + 482890".539) T 

+ 7". 455 T 2 + 0".008 T 3 (2.112) 


where Y = 360° = 1296000". 

With these fundamental arguments, the nutation quantities can then be represented by 


N ! 

■ 

* 1 1 

^=E 

(A 0 y + Ai jT) sin 

E k]t a i(T) I 

y=i 

■ 

4=1 ■* J 

N 


« 1 1 


(Boy + BiyT) cos 

E** a *( T ) 

y=i 


L »=i j j 


where the various values of a<, fcy,-, Ay , and By are tabulated in Table A. I. 


(2.113) 


(2.114) 


2.6. 2.1 Corrections to the 1980 IAU Model 


Additional terms can be optionally added to the nutations Sip and 6e in Eqs. (2.113) and (2.114). 
These include the out-of-phase nutations, the free-core nutations (Yoder, 1983) with period ujj (nom- 
inally 430 days), and the “nutation tweaks” A ip and As, which are arbitrary constant increments 
of the nutation angles Sip and Se. Unlike the usual nutation expressions, the tweaks have no time 
dependence. The out-of-phase nutations, which are not included in the IAU 1980 nutation series, are 
identical to Eqs. (2.113) and (2.114), with the replacements sin <-» cos: 


N 

- 

r s e 

sr = J2 

(A 2 j + A^jT) cos 


3 = 1 

. 



(2.115) 


N 

• 

r 5 y 

^° = E 

(i?2j + B^jT) sin 


3 = 1 

. 

i - 1 


(2.116) 
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Expressions similar to these are adopted fof the free-core nutations: 

Sip* = (Aoo + A 10 T) sin (w/T) + (A 20 A 30 T) cosfw/T) 

and 

Ss^ — {Boo + B l 0 T) cos(tL >/T) + (R 20 + B 30 T) sin(o;jT) (2.118) 

If the free-core nutation is to be retrograde, as expected on theoretical grounds, wj should be negative. 
The nutation model thus contains a total of 856 parameters: Aij (t=0,3; j = l,106) and Bij (i=0,3; 
j= 1,106) plus the free-nutation amplitudes i4 t o (i=0,3) , Bio (i=0,3). The only nonzero a priori 
amplitudes are the Ao/, Aiy, 5 0 y, Bij (j=l,106) given in Table A.I. 

The nutation tweaks are just constant additive factors to the angles Sip and 6 s: 


Sip —* Sip -f Aip 

(2.119) 

8 s — * 6 s + As 

( 2 . 120 ) 


Several alternatives are available as MODEST options to correct deficiencies in the IAU nutation 
model. The first possibility is to use empirically determined values of A ip t As as part of the polar 
motion and UT 1 input which was described in the next-to-last paragraph of section 2 . 6 . If this option 
is selected, the user is relying on nutation angles that are determined from other VLBI experiments 
near the date of interest, and performing linear interpolation. 

A second option employs the annual and semiannual amplitudes of Herring et al. (1986). These 
revised amplitudes axe given in Table VI in terms of the present notation, and in the units of Table 
A.I. 


Table VI 

Corrected Nutation Amplitudes (Herring et al. t 1986) 


Index, j 

9 

( 0 ". 0001 ) 

10 

( 0 ". 0001 ) 

Period, days 

182.6 

365.3 

In phase Ao,o 

-13172.2 

1471.0 

# 0,10 

5732.8 

72.1 

Out of phase ^ 2,9 

-8.3 

15.8 

- 82,10 

-2.9 

- 2.2 


Recent work by Zhu et al. (1989, 1990) has refined the 1980 IAU theory of nutation both by 
reexamining the underlying Earth model and by incorporating recent experimental results. The 
nutation series derived in that work are also available as MODEST modeling options. The Zhu et al. 
results are tabulated here in three parts: a) the original 106 terms of the 1980 IAU series with revised 
amplitudes in Table A. II, b) four sets of out of phase terms in Table A. Ill, and c) an additional 156 
terms due to planetary perturbations in Table A. IV. 

For simulation purposes, the older Woolard nutation model is also available in MODEST. With 
the exception of the number, amplitudes, and arguments of the terms, the older series is exactly 
analogous to the 1980 IAU theory, i.e., of the form of Eqs. (2.113) and (2.114). For completeness of 
documentation, the coefficients are listed in Table A.V. 

No partial derivatives with respect to Woolard or Zhu et al. amplitudes are currently calculated. 
It is emphasized that, for the present, the default nutation model in MODEST is just the 1980 IAU 
nutation model given in Table A.I. 
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2.6.3 PRECESSION 


The next transformation in going from the terrestrial frame to the celestial frame is the rotation 
jP. This is the precession transformation from mean equatorial coordinates of date to the equatorial 
coordinates of the reference epoch (e.gr., J2000). It is a composite of three rotations discussed in detail 
by Melbourne et al . (1968) and Lieske et al. (1977): 

( cos Z sin Z (A 

— sin Z cos Z 0 (2.121) 

0 0 Oj 


cos 0 0 sin 0 

Q(0) = ( 0 10 

— sin 0 0 cos 0 

cos £ sin f 0 
R(—$) = ( — sin £ cos £ 0 

0 0 1, 

P = R(-s)Q(e)R(-Z) 

cos £ cos 0 cos Z — sin £ sin Z cos £ cos 0 sin Z + sin £ cos Z cos £ sin 0 

— sin £ cos 0 cos Z — cos £ sin Z — sin f cos 0 sin Z + cos £ cos Z — sin f sin 0 

— sin 0 cos Z — sin 0 sin Z cos 0 

The auxiliary angles f , 0, Z depend on precession constants, obliquity, and time as 


( 2 . 122 ) 

(2.123) 

(2.124) 


f = 0".5mT + 0". 30188 T 2 + 0".017998 T 3 (2.125) 

Z = 0".5 mT + l". 09468 T 2 -f 0" .01 8 203 T 3 (2.126) 

0 = nT - 0". 42665 T 2 - 0".041833 T 3 (2.127) 


where the speeds of precession in right ascension and declination are, respectively, 

(2.128) 
(2.129) 


m = pls cos g 0 - PPL 
n = p /,5 sin £q 


and pls = the luni-solax precession constant, Ppl ~ planetary precession constant, £o = the obliquity 
at J2000, and T [Eq. (2.107)] is the time in centuries past J2000. Nominal values at J2000 are pls 
= 5038".7784/cy, ppl = 10".5526/cy; these yield the expressions given by Lieske et al. (1977) and 
Kaplan (1981): 

^ = 2306". 2181 T + 0". 30188 T 2 -f 0".017998 T 3 (2.130) 

0 = 2004". 3 109 T - 0". 42665 T 2 - 0".041833 T 3 (2.131) 

Z = 2306". 2181 T + l". 09468 T 2 + 0". 018203 T 3 (2.132) 


Partial derivatives of the VLBI observables with respect to luni-solar and planetary precession are 
derived from the expressions (2.124-2.129) and given in section 2.9. The precession matrix completes 
the standard model for the orientation of the Earth. Numerical checks of direct estimates of preces- 
sion corrections against similar estimates based on the perturbation rotation (next section) ensure 
consistency. 


2.6.4 PERTURBATION ROTATION 


This standard model for the rotation of the Earth as a whole may need a small incremental 
rotation about any one of the resulting axes. Define this perturbation rotation matrix as 


H = A x A y A* 


(2.133) 
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where 


1 

0 

0 


0 

50 , 

i 


(2.134) 


Ax 


0 

1 


-50, 


with 50, being a small angle rotation about the x axis, in the sense of carrying y into z; 


( i o -se v \ 

A y =10 1 0 (2.135) 

\SQ y 0 1 J 


with 50 y being a small angle rotation about the y axis, in the sense of carrying z into x; and 


A, = 




(2.136) 


with 50, being a small angle rotation about the z axis, in the sense of carrying x into y. For angles 
of the order of 1 arc second we can neglect terms of order 50 2 Rg as they give effects on the order of 
0.015 cm. Thus, in that approximation 


/ i 50 , -5e y ^ 

n = -50, 1 50, (2.137) 

V 5© tf —50* 1 J 

In general, 

6&i = 6Gi(t) = 6S i0 + SGiT + fi(T) (2.138) 

which is the sum of an offset, a time-linear rate, and some higher order or oscillatory terms. Currently, 
only the offset and linear rate are implemented. In particular, a non-zero value of 50 y is equivalent 
to a change in the precession constant. Setting 

50, =50 y = 50, =0 (2.139) 

gives the effect of applying only the standard rotation matrices. 

Starting with the Earth-fixed vector, ro, we have in sections 2.3 through 2.6 above shown how 
we obtain the same vector, r c , expressed in the celestial frame: 

r c = CIPNUXY(tq + A) (2.140) 


2.7 EARTH ORBITAL MOTION 

We now wish to transform these station locations from a geocentric celestial reference frame 
moving with the Earth to a celestial reference frame which is at rest relative to the center of mass 
of the Solar System. In this Solar System barycentric frame we will use these station locations to 
calculate the geometric delay (see Section 2.1). We will transform the time interval so obtained back 
to the frame in which the time delay is actually measured by the interferometer — the frame moving 
with the Earth. 

Let E' be a geocentric frame moving with vector velocity = fie relative to a frame, E, at rest 
relative to the Solar System center of mass. Further, let r(t) be the position of a point ( e.g ., station 
location) in space as a function of time, t , as measured in the E (Solar System barycentric) frame. 
In the E' (geocentric) frame, there is a corresponding position r'(t') as a function of time, t* . We 
normally observe and model r'(t') as shown in sections 2.3 through 2.6. However, in order to calculate 
the geometric delay in the Solar System barycentric frame (E), we will need the transformations of 
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r(t) and r'(t'), as well as of t and as we shift frames of reference. Measuring positions in units of 
light travel time, we have from Jackson (1975): 


00 

T'(S) = r(t) + {'i-l)r(t)ZZ-' 1 0t 

(2.141) 

t’ = l[t- r(t) 0} 

(2.142) 

and for the inverse transformation: 


r(t) = r'(t') + ( 7 -l)r'(t') ^ + 10t' 

(2.143) 


(2.144) 

where , 

(2.145) 


Let t\ represent the time measured in the Solar System barycentric frame (S), at which a wave 
front crosses antenna 1 at position i*i(£i). Let r 2 (£i) be the position of antenna 2 at this same time 
as measured in the Solar System barycentric frame. Also, let be the time measured in this frame 
at which that same wave front intersects station 2. This occurs at the position r 2 (t$). Following 
section 2.1, we can calculate the geometric delay t\ t%. Transforming this time interval back to the 
£' (geocentric) frame, we obtain 

*S # - = 7 ($ - h) - 7[r a(«5) - MM } * fi ( 2146 ) 

Assume further that the motion of station #2 is rectilinear over this time interval. That assumption 
is not strictly true but, as discussed below, the error made as a result of that assumption is much less 
than 1 cm in calculated delay. Thus, 



r 2 (t$) — r 2 (ti) + 02^2 ~ *0 

(2.147) 

which gives: 

r 2 (t 2 ) — r l(^l) = r 2(tl) — ri(tl) + 02^2 ~ *l) 

(2.148) 

and 

t *2 - t \ = ^(t 2 - ti) - ^[^(ti) - ri(tj)] 0 - 702 - 01*2 - *l] 

= "1(1 - 02 0){^2 — M “ 'y[ r 2( i l) - r l( t l)l ■ P 

(2.149) 


This is the expression for the geometric delay that would be observed in the geocentric (E') frame in 
terms of the geometric delay and station positions measured in the Solar System barycentric system 


Since our calculation starts with station locations given in the geocentric frame, it is convenient 
to obtain an expression for [r 2 (£i) — ri(^i)] in terms of quantities expressed in the geocentric frame. 
To obtain such an expression consider two events [r^ (t^), r 2 (^i )] that are geometrically separate, but 
simultaneous, in the geocentric frame, and occurring at time t\. These two events appear in the Solar 
System barycentric frame as: 

rx(ti) = r ;(ti) + (7 - l)ri(ti) ^ + 10t\ (2.150) 

and as: .. 

r 2 (i 2 ) = r'(t' x ) + (7 - + 1 fit\ (2-151) 

P 
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where 


(2.152) 


t2-h = l[r',(t’ 1 )-T' 1 (t' 1 )]fi 

With these three equations and the expression 

^ 2 (^ 2 ) = ^* 2 (^ 1 ) + $ 2(^2 ” ^l] (2.153) 


we may obtain the vector r 2 (ti): 

r 2 (ii) = r' (*') + { n - l)r' 2 (ti) • ^ - ^ 2 [r'(t'J -r'^)] 0 (2.154) 

P 

This is the position of station #2 at the time ti as observed in E. FYom this we obtain: 
ra(*i) -ri(tx) = r' ft) -r^) + ( 7 - l)(r'(^) - r^)] ■ ^ 

(2.155) 


As shown in section 2 . 1 , the vectors [r 2 (*i) — ri (ti)] and ft 2 aTe all * s needed to obtain t 2 — ti for 
the case of plane waves. For curved wave fronts we will need to know the individual station locations 
in the barycentric frame as well. These we obtain from (2.150) and (2.154) with t r x set equal to zero. 
Setting t\ — 0 is justified since the origin of time is arbitrary when we are trying to obtain time 
differences. 

In the actual coding of these transformations, the relationship for the transformation of velocities 
is also needed. Taking differentials of (2.143) and (2.144) we have: 


dr = dr' + (7 — l)dr' * 


00 

0 2 


+ ifidt' 


(2.156) 


dt = l{dt' + dr' ■ 0 ) 

Dividing to obtain dr/dt we obtain for station #2 in the E frame: 

02 + (7 ~ 1)02 -^ 2+70 

7 W+KH 


(2.157) 


(2.158) 


For station #2 relative to the geocentric origin, we have from (2.87) and ( 2 . 88 ): 


» (IPN 


dU 

dH 


XY r> E 


(2.159) 


where 


w E = 7.2921151467 x 10 -6 rad/sec 


(2.160) 


is the inertial rotation rate of the Earth as specified in Kaplan (1981), p. 12 . This is not a critical 
number since it is used only for station velocities, or to extrapolate Earth rotation forward for very 
small fractions of a day (i.e., typically less than 1000 seconds). Actually, this expression is a better 


dH 

approximation than it might seem from the form since the errors in the approximation, — - = 

dt 


are very nearly offset by the effect of ignoring the time dependence of PN. 

The assumption of rectilinear motion can be shown to result in negligible errors. Using the plane 
wave front approximation ( 2 . 2 ), we can estimate the error St in the calculated delay due to an error 
A ^ 2 m the above value of fi 2 : 


St = k [r 2 (fi) -ri(ti)] 


[l-k (0 2 + A0 2 ) 1 — k 02 J 


« rA0 2 


(2.161) 
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Further, from (2.158) above, 

A0 2 « A02 

since 

7 w 1+ 10“ 8 

For the vector £2 a f rame rotating with angular velocity u ;, 
time interval r due to neglecting the rotation of that frame is 

A#* » 


(2.162) 

(2.163) 

the error A/S^ that accumulates in the 


(2.164) 


Thus for typical Earth-fixed baselines, where r < 0.02 sec, neglect of the curvilinear motion of station 

#2 due to the rotation of the Earth causes an error of < 4x 10"“ 14 sec, or 0.0012 cm, in the calculation 

of r. Similarly, neglect of the orbital character of the Earth’s motion causes an error of the order of 
0.00024 cm maximum. 

The position, Re, and velocity, /3 Ei of the Earth’s center about the center of mass of the Solar 
System are: 

Re = (2.165) 

0E = - 5 ^- (2.166) 

where the index t indicates the Sun, Moon, and all nine Solar System planets, m, is the mass of the 
body indexed by t, while R * and fi i are that body’s center- of- mass position and velocity relative to the 
center of the Earth in the celestial frame. In a strict sense, the summation should be over all objects 
in the Solar System. Except for the Earth-Moon system, each planet mass represents not only that 
planet’s mass, but also that of all its satellites. The R* and are obtained from the JPL planetary 
ephemeris (DE200 as of May, 1982) for the J2000 frame. 

Working in a frame at rest with respect to the center of mass of the Solar System causes relativistic 
effects due to the motion of the Solar System in a “fixed frame” to be included in the mean position 
of the sources and in their proper motion. The effect of galactic rotation can be easily estimated. In 
the vicinity of the Sun, the period for galactic rotation is approximately 2.2 xlO 8 years. Our distance 
from the center is approximately 10 kpc = 3.086 xlO 22 cm. Thus, our velocity is 



27ri2 

T7 


« 9.3 x 10“ 4 


(2.167) 


For a source at zero galactic latitude, the maximum change in apparent position (over one half galactic 
rotation) is 


A0 « 


period 


« 5.5 x 10 6 arcsec/year 


(2.168) 


Since a 1-arcsecond angle subtends a distance of approximately 30 meters at one Earth radius, ne- 
glecting this effect is roughly equivalent to introducing an error of 0.015 cm/year on intercontinental 
baselines. For the present 12-year history of VLBI data, this implies a systematic error of the order 
of 0.2 cm. 
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2.8 ANTENNA GEOMETRY 


The above work indicates how the time delay model would be calculated for two points fixed 
with respect to the Earth’s crust. In practice, however, an antenna system does not behave as an 
Earth-fixed point. Not only are there instrumental delays in the system, but portions of the antenna 
move relative to the Earth. To the extent that instrumental delays are independent of the antenna 
orientation, they are indistinguishable to the interferometer from clock offsets and secular changes 
in these offsets. If necessary, these instrumental delays can be separated from clock properties by 
a careful calibration of each antenna system. That is a separate problem, treated as a calibration 
correction (e.g., Thomas, 1981), and will not be addressed here. 

However, the motions of the antennas relative to the Earth’s surface must be considered since 
they are part of the geometric model. A fairly general antenna pointing system is shown schematically 
in figure 5. The unit vector, s, to the apparent source position is shown. Usually, a symmetry axis 
AD will point parallel to s. The point A on the figure also represents the end view of an axis which 
allows rotation in the plane perpendicular to that axis. This axis is offset by some distance H from a 
second rotation axis BE. All points on this second rotation axis are fixed relative to the Earth. 
Consequently, any point along that axis is a candidate for the fiducial point which terminates this 
end of the baseline. The point we actually use is the point P. A plane containing axis A and perpen- 
dicular to BE intersects BE at the point P. This is somewhat an arbitrary choice, one of conceptual 
convenience. 



Figure 5. A generalized schematic representation of the geometry of a steerable antenna 


Consider the plane Q which is perpendicular to the antenna symmetry axis, AD, and contains 
the antenna rotation axis A. For plane wave fronts this is an isophase plane (it coincides with the 
wave front). For curved wave fronts this deviates from an isophase surface by « H 2 /[2R), where R 
is the distance to the source, and H is taken as a typical antenna offset AP. For H « 10 meters, 
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R = Rmoon — 60 Re » 3.6 x 10 8 meters, and the curvature correction H 2 j[2R) » 1.4 X 10~ 7 meters 
and is totally negligible. R has to be 5 km, or 10 _3 i? E , before this deviation approaches 1 cm 
contribution. Consequently, for all anticipated applications of radio interferometry using high-gain 
radio antennas, the curvature of the wave front may be neglected in obtaining the effect on the time 
delay of the antenna orientation. 

Provided the instrumental delay of the antenna system is independent of the antenna orientation, 
the recorded signal is at a constant phase delay, independent of antenna orientation, at any point on 
the Q plane. Since this delay is indistinguishable from a clock offset, it will be totally absorbed by 
that portion of our model. 
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Figure 0. Schematic representations of the four major antenna geometries used in VLBI 


2.8.1 Axis Offset 

The advantage of choosing the Q plane rather than some other plane parallel to it is that the axis 
A is contained in this plane, and the axis A is fixed relative to the BE axis by the antenna structure. If 
l is the length of a line from P perpendicular to the Q plane, the wave front will reach the Earth-fixed 
point P at a time A t = l/c after the wave front passes through the axis A. If r 0 is the model delay 
for a wave front to pass from P on antenna #1 to a similarly defined point on antenna #2, then the 
model for the observed delay should be amended as: 

r - to - (A t 2 - A*i) = t 0 + (/i - h)/c (2.169) 
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where the subscripts refer to antennas #1 and #2. 

For the inclusion of this effect in the model, we follow a treatment given by Wade (1970), Define 
a unit vector I along BE, in the sense of positive away from the Earth. Further, define a vector, L, 
from P to A. Without much loss of generality in this antenna system, we assume that s, L, and I are 
coplanar. Then: 


L = ±H 


fxjsxj) 

|lx[sxl] I 


(2.170) 


where the plus or minus sign is chosen to give L the direction from P to A. The plus sign is used if, 
when s and L are parallel or antiparallel, the antenna comes closer to the source as H increases. Since 


ix[sxi] = 8-1 [Isj (2.171) 

J = s-L = — [slf (2.172) 

where the sign choice above is carried through. 

Curvature is always a negligible effect in the determination of 8 • L. Likewise, gravitational 
effects are sufficiently constant over a dimension |L| so as to enable one to obtain to a very good 
approximation a single Cartesian frame over these dimensions. Consequently, it is somewhat easier 
to calculate a proper time At = l/c in the antenna frame and to include it in the model by adding it 
to ro, taking into account, in principle at least, the time dilation in going from the antenna frame to 
the frame in which r 0 is obtained. 


2*8.2 Refraction 


Thus, if So is the unit vector to the source from the antenna in a frame at rest with respect to the 
Solar System center of mass, perform a Lorentz transformation to obtain s, the apparent source unit 
vector in the Earth-fixed celestial frame. Actually, the antenna does not “look” at the apparent source 
position s, but rather at the position of the source after the ray path has been refracted by an angle e 
in the Earth’s atmosphere. This effect is already included in the tropospheric delay correction (Section 
4); however since the antenna model uses the antenna elevation angle Eq } the correction must be made 
here as well. For the worst case (elevation angle of 6°) at average DSN station altitudes, the deflection 
can be as large as 2 x 10“ 3 radians. Thus, 61 « He « 2 cm for H = 10 meters. A model option 
permits modification of sb to take atmospheric refraction into account. The large-elevation- angle 
approximation is the inverse tangent law: 

A E= 3.13 x 10“ 4 / tan E 0 (2.173) 

where E is the elevation angle, and A E the change in apparent elevation Eo induced by refraction. This 
model was implemented only for software comparison purposes, since it gives incorrect results at low 
elevation angles. In the notation of Section 4.2, a single homogeneous spherical layer approximation 
yields the bending correction in terms of the zenith troposphere delays pz , refractivity moment Mooi, 
scale height A, and Earth radius R: 


where 


AS = cos 1 [cos(S 0 + a 0 )/(l + Xo)] - «o 

(2.174) 

Xo = (PZi,„ + Pz wet /M 0 oi)/A 

(2.175) 

a 0 = cos -I [(l + a')/(l + a)] 

(2.176) 

<r = A/R 

(2.177) 

cr' — + <y[a + 2)/ sin 2 Eq^ — lj sin 2 Eq 

(2.178) 
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This formula agrees with ray- tracing results to within 1% at 6° and «15% at 1° elevation, while the 
corresponding comparisons for Eq. (2.173) give «25% at 6 and a factor of 3 at 1 . 

Since we are given I in terrestrial coordinates, we first perform the coordinate transformation 
given by Q above: 

l=Qiterre.trial 


With this done, obtain At = l/c , as shown in figure 6 for each of the major antenna types. Note 
that for “nearby” sources we also must include parallax ( e.g. f geographically separate antennas are 
not pointing in the same direction). If Ro is the position of the source as seen from the center of the 
Earth, and r is the position of a station in the same frame, then the position of the source relative to 
that station is 

R = Ro-r (2.180) 


and in (2.172) we make the substitution 



|Ro — r] • I 

l^o ~r\ 


(2.181) 


2.8.3 Unique Antennas 

One of the VLBI antennas employed by the IRIS project of the National Geodetic Survey does not 
fall into any standard category. It is unique because it is an equatorial mount designed for the latitude 
of Washington, D.C., but deployed at Richmond, Florida. The considerable latitude difference, and 
the axis offset of several meters, make it imperative that the antenna geometry be properly modeled. 
In the local VEN coordinate frame, the vector I is 

( sin 4>w A 

— cos 4>w € I (2.182) 

cos <f>w COS £ / 

Upon transformation to the Earth-fixed frame via the matrix VW [Eq. (2.67)], it becomes 

( cos A (sin <j>w cos — cos 4>w sin cos c ) + si n ^ cos 4>w si n € A 

sin X(s\n <f>w cos <j> — cos <j>w sin <f> cos c) — cos A cos <f>w sin e I (2.183) 

sin <f>w sin 4> + cos <f>y/ cos <j> cos c J 

Here (A, <ft) are the Richmond longitude and latitude, <f>w Is the latitude of Washington (39.06°), and 
e = 0.12° W of N is the azimuth misalignment. 

Two other one-of-a-kind antennas, Arecibo and Nancay, are seldom used in astrometric and 

geodetic VLBI work. The Arecibo antenna has hardware features which make it equivalent to an 

azimuth-elevation mount. The Nancay array has been treated by Ortega-Molina (1985), but the 
model is not presently incorporated in MODEST code. 

2.8.4 Site Vectors 

In the modeling software is the facility to provide a time-invariant offset vector in local geodetic 
coordinates (east, north, and local geodetic vertical) from this point (antenna location) to a point else- 
where, such as a benchmark on the ground. This is particularly useful in work involving transportable 
antennas which may be placed in slightly different places relative to an Earth-fixed benchmark each 
time a site is reoccupied. In modeling that offset vector, we make the assumption of a plane tangent 
to the geoid at the reference benchmark and assume that the local geodetic vertical for the antenna 
is parallel to that for the benchmark. With these assumptions there is an identity in the adjustments 
of antenna location with changes derived for the benchmark location. The error introduced by these 
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assumptions in a baseline adjustment is approximately A B X (d/i?#), where A B is the baseline ad- 
justment from its a priori value, d is the separation of the antenna from the benchmark, and Re is 
the radius of the Earth. To keep this error smaller than 0.01 cm for baseline adjustments of the order 
of 1 meter, d < 600 meters is required. 

More troublesome is that an error in obtaining the local vertical by an angle 60, when using 
an antenna whose intersection of axes is a distance, H , above the ground, can cause an error of 
HsinSQ « HSQ in measuring the baseline to the benchmark (Allen, 1982). Unless this error is 
already absorbed into the actual measurement of the offset vector, care must be taken in setting up 
the antenna so as to make 60 minimal. For a baseline error <0.1 cm, and an antenna height of 10 
meters, 60 < 20 arcseconds is required. Often plumb bobs are used to locate the antenna position 
relative to a mark on the ground. This mark is, in turn, surveyed to the benchmark. Even the 
difference in geodetic vertical from the vertical defined by the plumb bob may be as large as 1 arc 
minute, thus potentially causing an error of 0.3 cm for antennas of height 10 meters. Consequently, 
great care must be taken in these measurements, particularly if the site is to be repeatedly occupied 
by antennas of different sizes. 

2.8.5 Feed Rotation 

Another physical effect related to antenna structures is the differential feed rotation for circularly 
polarized receivers. Liewer (1985) has calculated the phase shift 9 for various antenna types. It is 
zero for equatorially mounted antennas. For altazimuth mounts, 

tan 6 = cos <f> sin h/(sin <f> cos 6 — cos <f> sin 6 cos h) (2.184) 

with <f> = station latitude, h = hour angle, and 6 = declination of the source. For X-Y mounts, two 
cases are distinguished: orientation N — S or E — W. The respective rotation angles are 

tan(— 9) — sin <f> sin h/{co&<f> cos 6 + sin ^sin 6 cos h) (N — S) (2.185) 

tan(-0) = — cos h/(sin 6 sin h) ( E — W) (2.186) 

The effect cancels for group delay data, but can be significant for phase delay data. The effect on 
phase delay is 

r=(*2-0i)// (2.187) 

where / is the observing frequency and 0* the phase rotation at station i. The feed rotation correction 
is now an optional part of the MODEST model. 

Finally, another small correction which accounts for the effect of orientation of HA-Dec and X-Y 
antennas on the tropospheric path delay was recently considered by Jacobs (1988). Details are given 
in the troposphere section, 4.4. 
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2.9 PARTIAL DERIVATIVES OF DELAY WITH RESPECT TO 
GEOMETRIC MODEL PARAMETERS 

With respect to any given parameter, the calculation of the time-delay model must be at least 
as accurate as the data is sensitive to that parameter. Consequently, such effects as the curvature of 
the wave fronts were considered. However, such detail is not necessary for determining the derivatives 
with respect to the relevant model parameters. Here, the plane wave approximation is sufficient. 
Iteration on the estimated parameters and the rapid convergence of an expansion of the time delay in 
the relevant parameters about some a priori point permit this simplification. 

In this plane wave approximation we wish to obtain the parameter derivatives with respect to: 

1. the nominal baseline components (actually, station locations), 

2 . the parameters of the whole Earth orientation matrix Q described in section 2 . 6 , 

3. the solid-Earth tidal parameters, 

4. the parameters of source location (right ascension and declination), 

5. the antenna axis offsets, 

6. the constant, 7 PPN) in the retardation of the light ray due to gravitational effects. 

The expressions for these derivatives are considerably simplified if tensor notation, with the 
Einstein summation convention, is employed. Before proceeding any further, we make the following 
definitions for this section: 

r = time delay modeled in the geocentric frame, 

t 9 = this same time delay, but modeled in the Solar System center of mass frame, 
b = source unit vector (in the celestial system at rest with respect to the Solar System 
center of mass), 

(} = velocity of the geocentric frame as measured in the Solar System center of mass frame 

(remember, all distances are measured in time; thus, this quantity is dimensionless), 

= velocity of station #2 in Solar System center of mass frame, 
p = 1 + s fi 2 - This is a factor » 1.0001, which arises from the motion of station #2 during 
the passage of the wave front from station #1 to #2, 

7 = (1-P ) f > 

72 = (I-P 2 ) > B t 

Q = matrix which transforms from the terrestrial system to the celestial system, 

Lo = the baseline vector in the terrestrial system, 

L, = this same baseline vector in the celestial system center of mass frame, 

L = this same baseline vector in the celestial system. 


With these definitions (2.149) may be written 



T = ^(l - ■ P 2 )t, ~ ifi ■ I*. 

(2.188) 

For plane waves from (2.2): 

k-[r 2 — ri ) M, _ t-L. 

1 — k fi 7 1 + 8 '^2 P 

(2.189) 

Thus, 

T — —7[1 — ~ 7 PkL> k 

(2.190) 

For parameters (represented 

symbolically by rj) associated with L $k only: 




(2.191) 

Define the vector: 

= -[7(1- /M.) ^f+ 7ft] 

(2.192) 
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Then 



dr) 


2.9.1 Source Parameters 

For parameters associated with the source position only: 


Since 


Define the vector: 
Then, 

For example: 
Then, 

and 

and 


drj p 


dsk _ Sk dp 
dr) p dr) 


P — 1 + 3 l^2i 




= --1(1 -A, ft)— 

P 


dsk _ Skfoj dsi 
dr) p dr) 


hl- 


*kP 2 t 


d$i 

dr) 


M ) ~ -7(1 “ft, Pi) 




6a- 


3 i02j 


ds 

da 


?L = m ?Zl 

dr, M] dr, 


s= [ cos 8 cos or, cos 6 sin ot } sin S ] 


= [ cos 6 sin Qt, cos 8 cos ot t 0 ] = [ A u A 2i ^3 ] 


ds 


— = [ — sin 6 cos a, — sin 6 sin a, cos 6 ] = [ Fi, F 2i F 3 ] 


Or, if we define the matrices: 

and 

then: 


dr 

5 - = MiAi 

da 

dr 

ds = M < F < 


Ax Fx 

G = | A 2 F 2 
A3 F3 

M = ( Mx,M 2 t M 3 ) 


dr dr 
da’ d6 


= MG 


(2.193) 


(2.194) 

(2.195) 

(2.196) 

(2.197) 

(2.198) 

(2.199) 

( 2 . 200 ) 
( 2 . 201 ) 

( 2 . 202 ) 

(2.203) 

(2.204) 

(2.205) 

(2.206) 
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For a linear model of source “proper motion” [Eqs. (2.85)-(2.86)], the partials of r with respect to 
the time rates of change of right ascension and declination (a, 5) are 

= (t- t 0 ) MG (2.207) 

where to is a reference time. 

2.9.2 Station Parameters 

For station location parameters the algebra is somewhat more complex. Since 

L, =r 3 (ti)-r 1 (ti) ( 2 - 208 ) 

= r 2 (t2) - ri(ti) - P 2 [h ~ <i] 

= r 2 (f2) ~ 1 ’ 1 (^ 1 ) — 'tfaMW) — r i (^i)l ' P 

= [r^(ii) - r + (7 - l)K(t'x) - ri(t'i)] ■ PP ~ 1 P2P • teW) ~ *!(*!)] 

we have: 

L, = L + (7 - 1)1 PP - 1P2P ' I* 

or in tensor notation _ 

L„ = *, + fi . L j 

Define the tensor: 

Eij = Si, + ^~ 2 1)A -1P2, Pi 

Then 

L g , = EijL } - 

Since 

Lj = QjkLo k 
L ti = EijQjkL o k 


(2.209) 

( 2 . 210 ) 

( 2 . 211 ) 

( 2 . 212 ) 

(2.213) 

(2.214) 


dr dr 
dd’ 38 


T = %EnQ jk L 0k (2.215) 

For parameters which are involved with station locations expressed in the terrestrial coordinate system: 

£-[***-1 < 2216 > 


where the vector element 

S fc = ViEijQi* (2.217) 

Such parameters are: r ap . (radius off spin axis), A? (longitude), z? (height above the equator), r api , 
X it z a (the station coordinates’ respective time rates), h 2i (vertical quadrupole Love number), hi 
(horizontal quadrupole Love number), 1 pi (phase lag of maximum tidal amplitude). The subscript 
refers to station number, z.e., i = 1, 2. Define the matrix: 


[— J?i, R 2 , — Ai, A 2 , — Z\ y Z 2 ,-R u R 29 — Ai, A 2 , — Zi 9 Z 2) —Vi, V 2 , —-Hi, H 2) ®i, $ 2 ] (2.218) 
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where each column contains the partials of the Lq component vectors x, y, 
parameters. For example, for the constant terms in the cylindrical station 
(2.38) through (2.40)]: 

( \ 


Ri = 


A, = 


dr %i 

d L 0y 

d_L o, 

V 3r? p , J 

( dLp . A 
aA° 



-r® p . sin A?' 
r »p. cos A< 


= 


dA? 
dLp, 

V 3a? y 

dL 0y 

dzf 
dLp, 

3*° y 

For the station coordinate rates, 

& = (t-t 0 )fl, Zi = (t 

From Eqs. (2.50) through (2.61), and relying on Williams (1970): 

/ \ 


Vi = 


2t 


#,= 


$i = 


dh 

d&iy 

dh 2 i 

d jii 

\ dh 2i J 

( dtix \ 

dhi 

dSiy 

dhi 

dSu 

V dhi j 

( dSix \ 
dip, 
dSi« 


= 5(i)V(,)W(.) 


= 5(t)y(t)W(.) 


t())Zi 



a 

d Sj, 

v. av< y 


= S(t)F(t)W(i) 


A) 

( ^gj 2) (*) A 

aVi 

5 ff2 2) (0 


av< 




(*) 


V J 


where % — 1 implies station #1, t = 2 implies #2, and 5(1) = —1, while 5(2) 
gW with respect to i/> are 


dip 


3 P,r% 

m 


hr p • R. 4 [t/pAT 9 XpV^] 


z with respect to the 
coordinates [see Eqs. 

(2.219) 

( 2 . 220 ) 

( 2 . 221 ) 

( 2 . 222 ) 

(2.223) 

(2.224) 

(2.225) 

= 1. These partials of 

(2.226) 
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dgffi _ 3 M« r p j l r pl 

d * R ‘ yfe + Vi 


[r p • R,][i p X + y P Y.\ - \x p Y. - y p X,] 2 


9gj, _ 3 **« r p 

d* Rf 


l \y P X, %pY t ] 


v/^p + j/ p -^« p -l2r p • R, z P 2,] 

* . / < t *2 _ i _ **2 


\j x l + i 


(2.227) 

(2.228) 


D = 


Also, define a vector: 

dr dr dr dr dr_ dr_ dr dr dr dr dr dr 
d^’ 9A°’ 9Ag’ «*i ’ 9 z 2’ ^Ai’ 3A 2 ’ 9*i ’ 9i 2 ’ 

3r dr dr dr dr dr dr dr 

dii ' dfc 2 * dh 2 i ’ dh .22 d/ 2 1 di 22 d^i dip? 


Then 


D = RVY 


(2.229) 

(2.230) 


2.9.3 Earth Orientation Parameters 


Certain parameters 
parameters 


such as UT1, polar motion, precession, 


— — • 



and nutation affect Q only. For these 


(2.231) 


Define a vector: 


= *kE ki 


(2.232) 


Then 



for parameters which affect only the orientation of the Earth as a whole. 


(2.233) 


2.9.3. 1 UT1 and Polar Motion 


A number of parameter partials are available for the orientation of the Earth. These are for 
UT1, X pole, Y pole, and nutation, as well as the angular offset and angular rate terms in the Earth 
orientation perturbation matrix Cl. From (2.87): 

Q = A PNUXY (2-234) 


Define the matrix: 



0 0 
0 — sin 02 

0 — cos ©2 


0 

COS ©2 I 

— sin ©2 j 


Then, the partial required for the Y polar motion parameter is: 


(2.235) 


- = n PNUXY 1 
9©y 


(2.236) 


An analogous technique is used for the X pole angle, working with the matrix partials . Partials 

with respect to UTl involve a slight complication due to the last three terms in Eq. (2.94). On the 
assumption that only the term linear in T u contributes significantly, 


dU 

d{UT 1) 


||(1 + 1/365.25) 


(2.237) 
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2.9.3. 2 Nutation 


Partial derivatives of the VLBI observables with respect to the nutation angles and amplitudes 
appear formidable at first sight, but are relatively easy to evaluate if the calculation is performed in 
an organized fashion. Symbolizing the parameters by 77, we need to evaluate the partials of the matrix 
Q with respect to 17: 


drj ! 


d6ip 


35$ / dr ) 


■XY 




(2.238) 

(2.239) 


dN 


Since Se — e — the first partial on the rhs of Eq. (2.239) is equal to . The derivatives of N 

de 

with respect to the angles Srp and Se are easily obtained from the expression for N in Eq. (2.105): 
dN 


dSrP 


— sin Srp cose cos sine cos Srp 

cosfcos^t^ — cos £ cos e sin Srp — cos £sin esin Srp 
— sin£cos£t/> — sinffcos esin£^ — sin Ssin esin Srp 


and 


dN 


— sin e sin Srp 


cos e sin Srp 


— — — I 0 — cos £ sin e cos 6 rp + sin t cos e cos e cos e cos Srp + sin £ sin e 

0 — sin £sin e cos Srp — cos £ cos e sin £ cos e cos Srp — cos £sine 

hVom Eq. (2.92), the partials of U with respect to Srp and Se are 


(2.240) 


(2.241) 


au 


dH 


dSrp ) Se 


— sin H — cob H 0 
cos H -sinff 0 
0 0 0 / dSlf, ’ Se 


and, from Eq. (2.97), 


dH 

— cose / (cos 2 Srp + cos 2 e sin 2 Srp) 
dH 

= — sine tan Srp / (1 + cos 2 e tan 2 Srp) 


(2.242) 

(2.243) 

(2.244) 


the {/-dependent terms in Eqs. (2.238) and (2.239) are evaluated. 

Partials of Srp and Se with respect to the parameters A{j and Bij are obtained immediately from 
Eqs. (2.117) and (2.118). For the “free nutations”, 


d6j>! 

3Aoo 

= sin ujf T } 

dSe ? 
dB 00 

= cos w/T 

(2.245) 

d6j>* 

dAiQ 

— T sin ojfT } 

dSe 1 

dB 10 

— T cos ojfT 

(2.246) 

dSrpf 

dA-zo 

= cos ujf T, 

d Se* 
dB2o 

— sin ojfT 

(2.247) 

d6il> 1 
< 5^30 

— T cos oj/Tj 

dSe? 

d#30 

= TsincjfT 

(2.248) 


and for the 1980 IAU series terms (j = 1 to 106): 
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(2.249) 


dSrp 

dAoj 


sin^^A: J iaj(T)j , 

t=i 


dSrp 

dAij 


5 

= rsin^fcy<aj(T)J , 

1=1 


dtnp 

dA 2 j 


t=i 


dSip 

dA zj 


— T cos [ y, fcyiQ, (T) j 


t=i 


dSe 

dB 0 j 

dSe 

dB u 

dje 

dBoj 

dSe 

dB^ 


= COS [y^,a,(T)] 

t= 1 

5 

= T cos [ y kjiOi (T)J 

i=i 

5 

= sin[y %a,(T)] 

i=l 

5 

= T sin Jy kjioti (T)j . 

i=i 


(2.250) 

(2.251) 

(2.252) 


2.9.3.S Precession 

Partial derivatives of the observables with respect to precession parameters are evaluated in a 
manner similar to those with respect to nutations. Symbolizing either the luni-solar precession p LS 
or the planetary precession ppl by tt, the partial of the rotation matrix Q is 

dP 

The partials are very complicated, and will be written in terms of the partials of each matrix 

d'K 

element Pa- 


— = — cos So sin f cos 0 cos Z/2 — sin So cos 5 sin 0 cos Z 

T dpLS 

— cos So cos £ cos 0 sin Z/2 — cos So cos f sin Z/2 

— cos So sin ( sin Z/2 (2.254) 


1 dPn 

T dppL 


~ sin $ cos © cos Z/2 4 cos f cos 0 sin Z/2 
4 cos £ sin Z/2 4 sin f cos Z/2 


1 &P 12 

T dpLS 


= — cos So sin f cos 0 sin Z/2 — sin So cos £ sin 0 sin Z 
4 cos Sq cos f cos 0 cos Z/ 2 + cos S 0 cos ( cos Z/2 


— cos Sq sin $ sin Z/2 


\ 3 P 

— = sin £ cos 0 sin Z/2 — cos £ cos 0 cos Z/2 

T 3 pp L 

— cos f cos Z/2 4 sin f sin Z/2 


1 dPi3 

T dpLS 


— cos So sin £ sin 0 /2 4 sin ?o cos £ cos 0 
^ - = Tsin f sin 0/2 


(2.255) 


(2.256) 


(2.257) 

(2.258) 
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1 

T dpLS 


— cos £ 0 cos f cos 0 cos Z/2 + sin £o sin £ sin © cos % 

+ cos So sin j cos 9 sin Z/ 2 4~ cos £o sin £ sin Z/2 

— cos Sq cos £ cos Z/2 


1 d P 

= cos £ cos 0 cos Z/2 — sin £ cos 0 sin Z/2 

r 

— sin £ sin Z/2 + cos f cos Z/2 


1 9J > 22 
T <5pl5 


cos £o cos f cos 0 sin Z/2 + sin £o sin f sin 0 sin Z 

cos So sin £ cos 0 cos Z/2 — cos £o sin f cos Z/2 
cos So cos £ sin Z/2 


1 5 JPOO * t m I 

— — = cos £ cos 0 sin Z/2 -f sin £ cos 0 cos Z/2 

T dp PL 

+ sin f cos Z/2 + cos f sin Z/2 


1 9^23 

r ^pls 


1 ^-^31 

r ^pl 5 


1 dP32 
T dpLS 


= — cos £0 cos f sin 0/2 — sin £0 sin ? cos 0 
= Tcos f sin 0/2 


<9^23 
dp pl 


= — sin £0 cos 0 cos Z -f cos £0 sin 0 sin Z/2 

= —T sin 0 sin Z/2 
dp pl ' 

= — sin to cos 0 sin Z — cos ffo sin © cos Z/2 

dPo 2 


dp pl 
d^33 
dpLS 


= T sin ©cos Z/2 
= —T sinffo sin© 


dP: 


33 


dppL 


= 0 


A check on the algebra may be performed by noting that 


3P 

dppL 





and 


dP dP dP 

= - cos to h T sin to —— 

dpLS dppi, do 

The corresponding partials of the U matrix are much simpler: 


dU 

dPLS 


I sin H 
— cos to I — cos H 

V 0 


cos H 0 \ 
sin H 0 j 

0 0 J 


(2.259) 

(2.260) 

(2.261) 

(2.262) 

(2.263) 

(2.264) 

(2.265) 

(2.266) 

(2.267) 

(2.268) 

(2.269) 

(2.270) 

(2.271) 

(2.272) 

(2.273) 
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du 

dppL 



/ COS So 


(2.274) 


2.9. 3.4 Rotational Tweaks 

Finally, the partials of the nutation matrix with respect to the ‘Hweaks” Arp and As are obtained 
by making the replacements (2.119) and (2.120) in N. and — are then seen to be identical 

to Eqs. (2.240) and (2.241), with the same replacements for Sip and Se. Expressions analogous to 
Eqs. (2.242) and (2.243) account for the shift of the equinox by nutation changes Sip and Se. If the 
o priori tweaks are zero, the partials are exactly equal to the expressions (2.240) and (2.241). 

For the parameters in the perturbation matrix, fl, from (2.138): 


an 

960*0 



0 o\ 
0 1 
-1 oj 


an 

96© z 



0 0\ 
0 t 
-to) 


(2.275) 


(2.276) 


where t is the number of years from the reference epoch ( t.g J2000). Then, by substituting these 
matrices for 0 in (2.138), we obtain the appropriate partials of Q for perturbations about the x axis. 
By analogy, the perturbation parameters about the y and z axes may also readily be obtained. 


2.9.4 Additive Parameters 

If we seek the partials of parameters that affect only the “add-on" terms in r = ro + Ar, then 
from (2.149) we have: 

( 2 - 277 ' 

for terms which were “added on” in the Solar System barycenter. An example is gravitational bending: 


For terms “added on” in the geocentric frame, then: 

dr dAr 
drj drj 

An example is the antenna offset vector. In this case: 

dr 

d (offset station #2) 


±\A-M 2 


(2.278) 


(2.279) 


(2.280) 


and 


dr 

^(offset station # 1) 


= ±\J i - [s i ] 2 


(2.281) 


where the choice of sign for each station is determined by the choice of sign for that station in the 
model portion. 
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SECTION 3 


CLOCK MODEL 

The frequency standards (“clocks”) at each of the two antennas are normally independent of each 
other. Attempts are made to synchronize them before an experiment by conventional synchronization 
techniques, but these techniques are accurate to only a few /xsec in epoch and » 10“ 12 in rate. 
More importantly, clocks often exhibit “jumps” and instabilities at a level that would greatly degrade 
interferometer accuracy. To account for these clock effects, an additional “delay” r c is included in the 
model delay, a delay that models the behavior of a station clock as a piecewise quadratic function of 
time throughout an observing session. Usually, however, we use only the linear portion of this model. 
For each station this clock model is given by: 

T C = Td + T c2 {t “ tref) + T c3 (t - t r ef) 2 /2 (3.1) 

The term, t re f y may be set by the user as a specific time (Julian date), or by default taken as the 
midpoint of the interval over which the a priori clock parameters, r cl , r c 2 , r c 3, apply. 

In addition to the effects of the lack of synchronization of clocks between stations, there are 
other differential instrumental effects which may contribute to the observed delay. In general, it is 
adequate to model these effects as if they were “clocklike”. However, the instrumental effects on delays 
measured using the multifrequency bandwidth synthesis technique (Thomas, 1981) may be different 
from the instrumental effects on delays obtained from phase measurements at a single frequency. 

The bandwidth synthesis process obtains group delay from the slope of phase versus frequency 

( d<f>\ 

T = across multiple frequency segments spanning the receiver passband. Thus, any instrumental 
contribution to the measured interferometer phase which is independent of frequency has no effect on 
the delay determined by the bandwidth synthesis technique. However, if delay is obtained directly 

from the phase measurement, <f> t at a given frequency, 1/, then this derived phase delay {^Tpd ” 
does have that instrumental contribution. V 

Because of this difference, it is necessary to augment the “clock” model for phase delay measure- 
ments: 

T c p d + ^"c4 [t ~ t re f ) + T c ${t — fre/) 2 / 2 (3-2) 

where r c is the clock model for bandwidth synthesis observations and is defined in (3.1). Since 
the present system measures both bandwidth synthesis delay and phase delay rate, all of the clock 
parameters described above must be used. However, in a “perfectly” calibrated interferometer, r c4 
= r c 5 = 0. This particular model implementation allows simultaneous use of delay rate data derived 
from phase delay with delay data derived by means of the bandwidth synthesis technique. However, 
our particular software implementation currently is inconsistent with the simultaneous use of delay 
derived from bandwidth synthesis and delay obtained from phase delay measurements. 

To model the interferometer delay on a given baseline, a difference of station clock terms is 
formed: 

^ ^ c 9ta.tion 2 ^‘ c »i4.tion 1 (3.3) 

Specification of a reference clock is unnecessary until the parameter adjustment step, and need not 
concern us in the description of the model. 

The partial derivatives of model delay with respect to the set of five parameters (r c i,T' C 2,7’ C 3,r C 4,r C 5) 
for each station are so trivial as to need no further explanation. 
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SECTION 4 


TROPOSPHERE MODEL 


In order to reach each antenna, the radio wave front must pass through the Earth’s atmosphere. 
This atmosphere is made up of two components: the neutral atmosphere and the ionosphere. In 
turn, the neutral atmosphere is composed of two major constituents: the dry and the wet. The dry 
portion, primarily oxygen and nitrogen, is very nearly in hydrostatic equilibrium, and its effects can be 
accurately estimated simply by measuring the barometric pressure. Typically, at sea level in the local 
zenith direction, the additional delay that the incoming signal experiences due to the troposphere is 
approximately 2 meters. Except for winds aloft, unusually strong lee waves behind mountains (e.g., 
Owens Valley, California), or very high pressure gradients, an azimuthally symmetric model based on 
measurements of surface barometric pressure is considered adequate. We have not yet investigated 
where this assumption breaks down, though “back-of-the-envelope” calculations indicate that, except 
in the unusual cases above, the error in such an assumption causes less than 1-cm error in the baseline. 

Unfortunately, the wet component of the atmosphere (both water vapor and condensed water in 
the form of clouds) is not so easily modeled. The experimental evidence (Resch, 1983) is that it is 
“clumpy”, and not azimuthally symmetric about the local vertical at a level which can cause many 
centimeters of error in a baseline measurement. Furthermore, because of incomplete mixing, surface 
measurements are inadequate in estimating this contribution which even at zenith can reach 20 to 
30 cm. Ideally, this tropospheric induced delay should be determined experimentally at each site. 
This is particularly true for short and intermediate (B < 1000 km) baselines, where the elevation 
angles of the two antennas are highly correlated in the observations. For long baselines, both the 
independence of the elevation angles at the two antennas, and the fact that often the mutual visibility 
requirements of VLBI constrain the antennas to look only in certain azimuthal sectors, allow the use 
of the interferometer data itself to estimate the effect of the water vapor as part of the parameter 
estimation process. For this reason, and because state-of-the-art water-vapor measurements are not 
always available, we also have the capability to model the neutral atmosphere at each station as a 
two-component effect, with each component being an azimuthally symmetric function of local geodetic 
elevation angle. 

At each station the delay experienced by the incoming signal due to the troposphere can be 
modeled using a spherical-shell troposphere consisting of a wet component and a dry component: 

Ttrop station < = T wet trop H" T dry trop ( 4 ’^) 


The total troposphere model for a given baseline is then: 


— Tt r op station 2 ^trop station 1 

If Ei is the apparent geodetic elevation angle of the observed source at station t, we have (dropping 
the subscript t): 

Ttrop = PZ drv Rdry{E) + PZ wet Rwet (#) ( 4 * 3 ) 

where pz is the additional delay at local zenith due to the presence of the troposphere, and R is an 
elevation angle mapping function. 

For some geodetic experiments, the observed delay has been corrected for the total tropospheric 
delays at the two stations, which are in turn calculated on the basis of surface pressure measurements 
for the dry component, and water- vapor radiometer measurements for the wet component. This 
correction is recorded in the input data stream in such a way that it can be replaced by a new model. 
In the absence of such external calibrations, it was found that modeling the zenith delay as a linear 
function of time improves troposphere modeling considerably. The dry and wet zenith parameters are 
written as 

PZi,* = P°Zi, m + PZi.vi* - *o) ( 4 ‘ 4 ) 

where to is a reference time. 
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Since the model is linear in the parameters p° and p, the partial derivatives with respect to zenith 
delays and rates are trivial. They are: 


dr 


d p%< 




(4.5) 


and 


dr 


dpz u 


or id 




where /(»’) = 1 for station #2, and —1 for station #1. 


4.1 CHAO MAPPING FUNCTION 


( 4 . 6 ) 


The simplest mapping function implemented in MODEST code is that obtained by C. C. Chao 
(1974) through analytic fits to ray tracing, a function which he claims is accurate to the level of 1% 
at 6° elevation angle and becomes much more accurate at higher elevation angles. 


where 


R = 


sin 2? + 


tan E + B 


A dry = 0.00143 
B d ry = 0.0445 
A wct = 0.00035 
B wet = 0.017 

The user must specify values for the zenith delays. 

The partial derivatives of delay with respect to the parameters A dry and B dry are: 


dr 

dA dry 


-/(*>Zirv R dry ! (tan E + Bdry) 


(4.7) 


(4.8) 

(4.9) 

(4.10) 

(4.11) 


(4.12) 


and 


dr 

&B dry 


f{')pZ dr yR d ryAdrv/{tznE+ B dr y ) 2 


where R dr y is the Chao mapping function, and E is the elevation angle. 


4.2 LANYI MAPPING FUNCTION 


(4.13) 


Analyses of intercontinental data indicate that the Chao mapping function [Eq. (4.7)] is inade- 
quate. To rectify this situation, two modifications have been made to the MODEST code. First, the 
dry-troposphere mapping parameters A dry and B dry of the Chao mapping function R dry have been 
promoted to the status of estimable parameters. Second, the code now permits the use of two more 
accurate mapping functions. The first of these is the analytic function developed by Lanyi ( 1984 ). In 
its simplest form, this mapping function employs average values of atmospheric constants. Provision is 
made for specifying surface meteorological data acquired at the time of the VLBI experiments, which 
may override the average values. Using numerical fits to ray-tracing results, Davis et al. ( 1985 ) have 
arrived at another function, designated the CfA-2.2 mapping function. Comparisons indicate that the 
Lanyi and CfA functions are in agreement to better than 1 cm over an extreme range of atmospheric 
conditions down to 6° elevation angles. Finally, an approximate partial derivative is obtained with 
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respect to one parameter in the Lanyi mapping function; this permits adjustment even in the absence 
of surface data. The Lanyi function was made the default MODEST troposphere model in early 1986. 

Motivation for and full details of the development of a new tropospheric mapping function are 
given by Lanyi (1984). Here we attempt to give a minimal summary of the final formulas. The 
tropospheric delay is written as: 

ftTop — F(E)/ sin E (4-14) 


where 

F[E) = pz dry F dry {E) + pz W ',F we t(E) 

+ \p\ iry W) + 2 p Zi , yP z„ME) + p 2 Zw ' t F b z(E))/A + p% irv F bi (E)/A 2 (4.15) 

The quantities p ZiTJI and p Zwet have the usual meaning: zenith dry and wet tropospheric delays. A 
is the atmospheric scale height, A = kT 0 /mg c , with k = Boltzmann’s constant, T 0 = average surface 
temperature, m = mean molecular mass of dry air, and g c = gravitational acceleration at the center o 
gravity of the air column. With the standard values k = 1.38066 X 10“ 16 erg/K, m = 4.8097 x 10 
g, gc — 978.37 cm/sec 2 , and the average temperature for DSN stations To = 292 K, the scale height 

A = 8567 m. ( . 

The dry, wet, and bending contributions to the delay, Fdry{F) } F we t(E)> and Fbi,b2,b3M\F)j are 


expressed in terms of moments of the refractivity as 

Fdry{E) = i4io (£) G ( AMno , u ) + 3 < 7uM 2 io <? 3 ( Muo, u )/ 2 ( 416 ) 

F wet (E) = A 0 i ( S ) G ( AM 10 i / Mooi , «)/ M 00 i ( 4 - 17 ) 

Fbi{E) = [crG 3 (JWuo, u)/sin 2 E - M 0 2oG 3 (Mi2o/Afo20, u )]/2 tan 2 E ( 4 - 18 ) 

F b 2 {E) = -MoiiG 3 (M n i/Moii,u)/2Mooitan 2 ^ ( 4 -19) 

F b a{E) — -A/oo 2 G 3 (Afio 2 /M) 02 , «)/2A^ooi tan ^ ( 4-2 ®) 

F U {E)= Afo3oG 3 (M 130 /Mo3o» u )/ tan 4 E ( 4 -2l) 

A misprinted sign in the last of Eqs. (5) of Appendix B of Lanyi (1984) has been corrected in Eq. 
(4.21). Here G(q,u) is a geometric factor given by 

G{q, u) = (1 + gu) _1/2 ( 4 - 22 ) 


with 


u = 2<j j tan 2 E (4.23) 

where a = A/R is a measure of the curvature of the Earth’s surface with standard value 0.001345. 

The quantities A, m (£) and M ilm are related to moments of the atmospheric refractivity, and are 
defined below. A l0 {E) involves the dry refractivity, while A 0 i{E) is the corresponding wet quantity. 
The Ai m {E) are given by 


Aim (-E) — -Woim + 2"Jt!(n— Jfc)! 

n= 1 Jt = 0 V 7 



U 

n 

\Mum 


1 + XuMum/Moim 


Mo lm 


(4.24) 


with the scale factor A = 3 for E < 10° and A = 1 for E > 10°. Only the two combinations (/, m) = 
(0,1) and (1,0) are needed for the A lm (E). The moments of the dry and wet refractivities are defined 


as 


oo 

M ni j = J dq q n f' dry {q)fL t {g) 


(4.25) 


where f dry , wet (q) are the surface-normalized refractivities. Here, n ranges from 0 to 1, i from 0 to 3, 
and j from 0 to 2; not all combinations are needed. Carrying out the integration in Eq. (4.25) for a 
three-section temperature profile gives an expression for the general moment 


M nij /n\ = (1 


■ e~ aqi )/a n+1 + e - ° 91 


[l-T 2 fc+n+1 (qi,92)]n 

L 1=0 


6 + i + l 


+ e -a, i T 6+n+ 1 ( g i, (?2 )/ a n+ 1 


(4.26) 
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Here, 


22 (^ 1 , 92 ) = 1- (92 ~9i)/« 


(4.27) 


The quantities qx and q 2 are the scale-height normalized inversion and tropopause altitudes, respec- 
tively. For the standard atmospheric model, qx = 0.1459 and g 2 = 1.424. The constants a and b are 
functions of the dry (a = 5.0) and wet (/? = 3.5) model parameters, as well as of the powers of the 
refractivities (» and j) in the moment definitions. Table VII gives the necessary a’s and 6’s. 


Table VII 

Dependence of the Constants a and 6 
on Tropospheric Model Parameters 


t 

3 

a 

6 

1 

0 

1 

a — 1 

0 

1 

p 

a/0-2 

2 

0 

2 

2(a-l) 

1 

1 

P+l 

0{ot + 1) - 3 

0 

2 

20 

1 

cT 

3 

0 

3 

3(a- 1) 


Note that the normalization is such that A /010 = this moment has therefore not been explicitly 
written in Eqs. (4.16) through (4.21). 

At present, provision is made for input of four meteorological parameters to override the default 
(average) values of the Lanyi model. These are: 1) the surface temperature To, which determines the 
atmosphere scale height (default value 292 K); 2) the temperature lapse rate IV, which determines 
the dry model parameter a (default values W = 6.8165 K/km, a = 5.0); 3) the inversion altitude 
hi, which determines q\ — hi/A (default value h\ = 1.25 km); and 4) the tropopause altitude / 12 , 
which determines <72 = ^/A (default value h <2 = 12.2 km). A fifth parameter, the surface pressure 
po, is not used at present. Approximate sensitivity of the tropospheric delay (at 6° elevation) to the 
meteorological parameters is —0,7 cm/K for surface temperature, 2 cm/(K/km) for lapse rate, and 
—2 cm/km for inversion and 0.5 cm/km for tropopause altitude, respectively. 

Partials of the delay with respect to the dry and wet zenith delays are obtained from Eqs. (4.14) 
and (4.15): 


— f[i)[F dr y{E) + 2pz drv Fbi{E)/A]/ sin E 
+ [ 2 Pz„. t Fb2(E)/& + 3p 2 Ziry F U {E)/ A 2 ]/ sin E (4.28) 

= f{i)[Evact[E) + 2pz iry Fb2(E)/A +2pz y ,«Fb3{E)/A]/ sinE (4.29) 


In analysis of data for which meteorological parameters are not available, it is convenient to introduce 
an approximation to the mapping function [Eqs. (4.14) and (4.15)] which involves a one-parameter 
estimate. This parameter p accounts for deviations from standard meteorological conditions. The 
tropospheric delay is expressed as 


Ttrop = {pz dTy + Pz wet )/ sin E + p 


dp 


(4.30) 
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where the partial derivative is 

dTtrop [PZj +PZ«,«t)»-^110 

dp G[Mno, w)[l + <7(Mno, «)] sin E 

pz„„_ t u [Mno - Aiioi/Mpoi) 

+ G(M no> «)G(Mioi/M 00 i, «) [G(M no , «) + G(M 10 i/Mooi, «)] sin E 


(4.31) 


4.3 CfA MAPPING FUNCTION 


Another approach to improved modeling of tropospheric delay was published by Davis et al. 
(1985). Analytic fits to ray- tracing results yield the CfA-2.2 mapping function 


R = 


sin E + 


tan E + — 


sin E + c 


(4.32) 


where E is the elevation angle. The three parameters a, 6, c are expressed in terms of meteorological 
data as 


a = 0.0002723 [ 1 + 2.642 X 10 _4 p 0 - 6-400 x 10 -4 e 0 + 1-337 x 10 -2 T o 

- 8.550 X 10 _2 a - 2.456 X I0~ 2 h 2 ] (4.33) 

b = 0.0004703 [ 1 + 2.832 x 10 -5 p 0 + 6.799 x 10 _4 « 0 + 7.563 x 10 _3 T o 

- 7.390 X 10 _2 a - 2.961 X 10 _2 h 2 ] (4-34) 

c = - 0.0090 (4-3 5 ) 


Here, po is the surface pressure and «o the surface partial water vapor pressure, both measured in 
millibars. The quantities T 0 , a, and h 2 have the same meaning and units as in Section 4.2. This 
function is one of three optional mapping functions in the MODEST model. In connection with 
testing parameter estimation for the Lanyi function, the partial derivative of delay with respect to 
surface temperature To in the CfA-2.2 function was also evaluated. It is 

d T _ PZjr^ry [ 3 - 64 l x 10~ 6 (sinE+ c)[tan + 6/(sin E + c)] — 3.557 X 10~ 6 a] ^ 

dTo (sin E + c) [tan E + 6/(sin E + c)\ 2 


4.4 ANTENNA AXIS OFFSET ALTITUDE CORRECTION 

Antennas with non-zero axis offsets, whose second rotation axis (A in figure 5) moves vertically 
with changing orientation, have zenith troposphere delays that may vary by 1 to 2 mm. Equatorial and 
X-Y mounts fall in this class (see figure 6). At low elevation angles this zenith variation is magnified by 
the mapping function to 1-2 cm. These variations must be modeled in experiments whose accuracies 
are at the millimeter level (e.g. short-baseline phase delay measurements). Memoranda by Jacobs 
(1988, 1991) derive the corrections based on considering only the dry troposphere component, and 
including all terms necessary to achieve an accuracy of a few millimeters. The correction to be added 
to the zenith dry tropospheric delay is 

Sr = -pz t ' V (H/A) fP (4-37) 

where H is the antenna axis offset, A the dry troposphere scale height (» 8.6 km), and ip is an angular 
factor that varies with the type of mount. For equatorial mounts, 

= cos <f> cos h (4.38) 
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where <p is the geodetic latitude and h the local hour angle east of the meridian. The Richmond 
antenna correction has this form with <f> replaced by <f>w and h by a pseudo-hour angle h R (see Section 
2.8.3), where 

h R — arctan |cos i?sin(0 — e)/[cos <pw sin E — sin <f>w cos Ecos(0 + c)]j (4.39) 

For north-south oriented X-Y mounts, 

xp = Bin E/(l — cos 2 & cos 2 E) 1 / 2 (4.40) 

where E is the elevation angle and 6 the azimuth (E of N). Finally, for east- west oriented X-Y mounts, 

xp = sin£'/(l — sin 2 6 cos 2 E) 1 ^ 2 (4.41) 
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SECTION 5 


IONOSPHERE MODEL 


The second component of the Earth’s atmosphere, the ionosphere, is a layer of plasma at about 
350 km altitude, created primarily by the ultraviolet portion of the sunlight. In the quasi-longitudinal 
approximation (Spitzer, 1962) the refractive index of this medium is 




— I! 1/2 


where the plasma frequency, v p , is 


v P = (pc 2 r 0 /7r) 1/2 « 8.97 X 10 3 p 1/2 


the electron gyrofrequency, v g) is 


Vn = 


eB 

2irmc 


(5.1) 


(5.2) 


(5.3) 


and © is the angle between the magnetic field B and the direction of propagation of the wave front. 
Here p is the number density of the electrons, and ro is the classical electron radius. 

For the Earth’s ionosphere, with p « 10 12 electrons/m 3 , i/ p « 8.9 MHz, while for the interplane- 
tary medium with p » 10 7 - 10 8 electrons/m 3 , v p w 28 - 89 kHz. In the interstellar medium, p w 10 
electrons/m 3 , which gives v v m 3 kHz. At typical microwave frequencies used for geodetic VLBI (8.4 
GHz), v p ju = 10- 3 for the ionosphere, 10“ 5 for the interplanetary medium, and 3 x 10 for the 
interstellar medium. 

The gyrofrequency, v gy for an electron in the « 0.2 gauss field of the Earth is * 0.6 MHz, Thus, 
for the ionosphere, vjv « 2 x HT 4 at S band (2.3 GHz), and uju « 7 X HT 5 at X band (8.4 GHz). 
For the interstellar medium B » 1CT 6 gauss, while for the interplanetary region B « 10 gauss. 

Relative to vacuum as a reference, the phase delay of a monochromatic signal transiting this 
medium of refractive index n is 


Tpd 


= ;/(- -e / Cv) l 1 + 1 (t)' + Kv)’ + -1* ,5 ' 4) 


where 


,- 1/2 


For 8.4 GHz, we may approximate this effect to parts in 10 6 — 10 7 by: 


A 


t-l 


1± 


(?) 


cos 0 




-q 


-e) 


cos 0 


(5.5) 


(5.6) 


where 




cro/« 

2x 


(5.7) 


and where I e is the t otal numbe r of electrons per unit area along the integrated line of sight. If we 
also neglect the term {v g cos®)/v, then the expression for A p<f becomes simple and independent of 
the geometry of the traversal of the wave front through the ionosphere: 


A P d = -q/v 2 


(5.8) 


This delay is negative. Thus, a phase advance actually occurs for a monochromatic signal. Since phase 
delay is obtained at a single frequency, observables derived from phase delay ( e.g ., phase delay rates) 
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experience an increment which is negative (the observable with the medium present is smaller than it 
would be without the medium). In contrast, group delays measured by a technique such as bandwidth 


f d<p v 

synthesis (r = — ) experience an additive delay which can be derived from (5.8) by differentiating 
<j> = vApd with respect to frequency: 


&gd = qfv 2 


(5.9) 


Notice that the sign is now positive, though the group delay is of the same magnitude as the phase 
delay advance. For group delay measurements, the measured delay is larger with the medium present 
than without the medium. 

For a typical ionosphere, r « 1 - 20 X 10“ 10 sec at local zenith for v = 8.4 GHz. This effect has 
a maximum at approximately 1400 hours local time and a broad minimum during local night. For 
long baselines, the effects at each station are quite different. Thus, the differential effect may be of 
the same order as the maximum. 

For the interplanetary medium and at an observing frequency of 8.4 GHz, a single ray path 
experiences a delay of approximately 6 x 10“ 7 sec in transiting the Solar System. However, the 
differential between the ray paths to the two stations on the Earth is considerably less, since the 
gradient between the two ray paths should also be inversely proportional to the dimensions of the 
plasma region (e.g., one astronomical unit as opposed to a few thousand kilometers). The ray path 
from a source at a distance of 1 megaparsec (3 X 10 7 km) experiences an integrated plasma delay of 
approximately 5000 seconds for a frequency of 8.4 GHz. In this case, however, the typical dimension 
is also that much greater, and so the differential effect on two ray paths separated by one Earth radius 
is still not as great as the differential delays caused by the Earth’s ionosphere. 


5.1 DUAL-FREQUENCY CALIBRATION 


These plasma effects can best be removed by the technique of observing the sources at two 
frequencies, i/i and v 2l where and where \u 2 - ^i |/(^2 + i'i) « 1. Then at the two 

frequencies v i and v 2 we obtain 

Tvl = T + qfi/l (5.10) 

and 

r„ 3 = t + q/v\ (5.11) 

Multiplying each expression by the square of the frequency involved and subtracting, we obtain 


t = ar u 2 + 


(5.12) 


where 


"f - v \ 


and 


b = 


A - A 


(5.13) 


(5.14) 


This linear combination of the observables at two frequencies thus removes the charged particle con- 
tribution to the delay. 

For uncorrelated errors in the frequency windows, the overall error in the derived delay can be 
modeled as 

^ = + (5.15) 


Modeling of other error types is more difficult and will not be treated in this report. Since the values 
of a and b are independent of q , these same coefficients apply both to group delay and to phase delay. 

If we had not neglected the effect of the electron gyrofrequency in the ionosphere, then instead 
of (5.12) above, we would have obtained 


r = &t U 2 + br^i *f 


q i/ g cos 0 

^ 2^1 (^2 - ^ l ) 


( 5 . 16 ) 


56 


where a and b are defined as in (5.13) and (5.14), respectively. 

If we express the third term on the right-hand side in units of the contribution of the ionosphere 
at frequency v 2 , we obtain 


t = ar U 2 + 1>t v 1 + 


Ap<jI /2 V g COS 0 

Vlfo + ^l) 


(5.17) 


For X band A pd « 1 - 20 xlO -10 sec at the zenith. When using S band as the other frequency in the 
pair, this third term is « 2 x 10~ 4 A pd cos@ « 2 - 40 xlCT 13 sec at zenith. In the worst case of high 
ionospheric electron content, and at low elevation angles, this effect could reach 0.1 cm of total error 
in determining the total delay using the simple formula (5.12) above. Notice that the effect becomes 
much more significant at lower frequencies. 

In the software chain used at JPL, the dual-frequency correction is performed prior to the process- 
ing step “MODEST” (Lowe, 1991). MODEST does not have the facility to perform this correction. 
However, the process is described here because it is important to understanding the data input to 
MODEST. For millimeter accuracy, or for lower observing frequencies even at centimeter accuracy 
levels, a correction for the gyrofrequency effect is necessary. 


5*2 TOTAL ELECTRON CONTENT 


In the absence of the dual-frequency observation capability described above, one can improve the 
model of the interferometer by modeling the ionosphere, using whatever measurements of the total 
electron content are available. The model we have chosen to implement is very simple. Its formalism 
is very similar to that of the troposphere model, except that the ionosphere is modeled as a spherical 
shell for which the bottom is at the height hi, above the geodetic surface of the Earth, and the top of 
the shell is at the height / 12 , above that same surface (see figure 7). For each station the ionospheric 
delay is modeled as 

n=kgI e S{E)/v 2 (5.18) 


where 


. 0.1cr o 

k = — 

2tt 


(5.19) 


I e is the total electron content at zenith (in electrons per meter squared xl0~ 17 ), and g = l(— 1) for 
group (phase) delay. E is the apparent geodetic elevation angle of the source, S{E) is a slant range 
factor discussed below, and v is the observing frequency in gigahertz. 

The slant range factor (see figure 7) is 


JR 2 sin 2 E + 2 Rh 2 + t% - Jr 2 sin 2 E + 2 Rh x + h\ 

S{E) = 7 f — (5.20) 

ft 2 — rtl 

This expression is strictly correct for a spherical Earth of radius R y and a source at apparent elevation 
angle E. The model employed uses this expression and a geoid surface with a local radius of curvature 
at the receiving station of R equal to the distance from the receiving station to the center of the Earth. 
The model also assumes this same value of R can be used at the ionospheric penetration points, e,g.: 


Ri = R + hi 


(5.21) 


This is not strictly true, but is a very close approximation, particularly compared to the crude nature 
of the total electron content determinations on which the model also depends. The total ionospheric 
contribution on a given baseline is 


7/ Citation 3 Ti station 1 (5.22) 

We assume that the ionospheric total electron content, J e , is the sum of two parts, one obtained by 
some external set of measurements such as Faraday rotation or GPS techniques, and the other by 
some specified additive constant: 

Ie = Ie meaa + h add (5.23) 
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These external measurements, in general, are not along directions in the ionosphere coincident with 
the ray paths to the interferometer. Thus, for each antenna, it is necessary to map a measurement 
made along one ray path to the ray paths used by the interferometer. Many different techniques to 
do this mapping have been suggested and tried; all of them of dubious accuracy. In the light of these 
problems, and in the anticipation that dual-frequency observations will be employed for the most 
accurate interferometric work, we have implemented only a simple hour-angle mapping of the time 
history of the measurements of I c at a given latitude and longitude to the point of interest. In this 
model we allow the user to adjust the “height”, h, of the ionosphere, but require 

hi — h — 35 km 

h 2 = h + 70 km (5.24) 



Figure 7. The geometry of the spherical ionospheric shell used for ionospheric corrections 

Nominally, this “height” is taken to be 350 km. Setting this height to zero causes the program to 
ignore the ionosphere model, as is required if dual- frequency observations have already been used to 
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remove the plasma effects. As in the troposphere model, these corrections can also be incorporated 
into the input data stream. Then the user is free to accept the passed correction, and use this model 
as a small alteration of the previously invoked model, or to remove the passed corrections. 

The deficiencies of these ionosphere models for single-frequency observations are compounded 
by the lens effect of the solar plasma. In effect, the Solar System is a spherical plasma lens which 
will cause the apparent positions of the radio sources to be shifted from their actual positions by an 
amount which depends on the solar weather and on the Sun-Earth-source angle. Since both the solar 
weather and the Sun-Earth-source angle change throughout the year, very accurate observations over 
the time scale of a year will be virtually impossible. 

Only one parameter is present in the ionosphere portion of the model. Again, the model is linear 
in the parameter I e add- Thus, the partial derivative with respect to this parameter is 

dr _ k /(station #) gr(data type) S[E) ^ 25) 

die add 1/2 

with /( 2) = 1 and /(l) = —1. 
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SECTION 0 


MODELING THE PHASE DELAY RATE (FRINGE FREQUENCY) 


The interferometer is capable of producing several data types: group delay, phase delay, and the 
time rate of change of phase delay. Actually, the time rate of change of group delay is also available. 
However, it is not accurate enough to be of significance for geodetic uses. The models discussed above 
are directly applicable either to group delay or to phase delay. However, the model for the time rate 
of change of phase delay (fringe frequency) must be either constructed separately, or its equivalent 
information content obtained by forming the time difference of two phase delay values constructed 
from the delay-rate measurements as shown below. We chose the latter route since then only models of 
delay are needed. The two phase delay values, r pc i(t±A), used to represent the delay-rate measurement 
information content are obtained from the expression 

Tpd{t =t A) = r m (t i A) + T r (t) i T r A (6*1) 

where r m (t) is the model used in the delay extraction processing step, r r {t) is the residual of the 
observations from that model, and f r is the residual delay rate of the data relative to that model. 
This modeling for the delay extraction step is covered in Thomas (1981), and is done in analysis 
steps prior to and completely separate from the modeling described in this report. The output of 
those previous steps is such that the details of all processing prior to the modeling described here 
are transparent to this step. Only total interferometer delays and differenced total interferometer 
phase delays (these phase delays are divided by the time interval of the difference) are reported to 
this step. One of the requirements of these previous processing steps is that the model delay used 
be accurate enough to provide a residual phase that is a linear function of time over the observation 
interval required to obtain the delay information. A linear fit to this residual phase yields the value 
of T r , the residual delay rate. Using these two values of T pd} obtained by (6.1) above, the quantity, R } 
is constructed by the following algorithm: 

R = [ r p<*fo + A) — r p( j[t — A)] 

2A t 6 ' 2 ' 

This value and the group delay measurement, r gdi are the two data types that normally serve as 
the interferometer data input to be explained by the model described in this report. The software, 
however, also has the option to model phase delay, r pd , directly. In the limit A — ► 0, this expression 
for differenced phase delay approaches the instantaneous time rate of change of phase delay (fringe 
frequency) at time t. In practice, A must be large enough to avoid roundoff errors that arise from 
taking small differences of large numbers, but should also be small enough to allow R to be a reasonably 
close approximation to the instantaneous delay rate. A suitable compromise appears to be A « 2 
seconds. Fortunately, A has a wide range of allowed values, and the capability to model interferometer 
performance accurately is relatively insensitive to this choice. 
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SECTION 7 


PHYSICAL CONSTANTS USED 


In the software that has been implemented we have tried to use the constants recommended by 
the IAU project MERIT (Melbourne et a/., 1983). Those that have not been defined in the text above, 
but which have an effect on the results that are obtained using the JPL software, are given below: 


Symbol 

Value 

Quantity 

c 

299792.458 

Velocity of light (km/sec) 

ro 

2.817938 x 1CT 15 

Classical radius of the electron (meters) 

Re 

6378.140 

Equatorial radius of the Earth (km) 


7.2921151467 x 10“ 5 

Rotation rate of the Earth (rad/sec) 

f 

298.257 

Flattening factor of the geoid 

h 2 

0.609 

Vertical quadrupole Love number 

1 2 

0.0852 

Horizontal quadrupole Love number 

hz 

0.292 

Vertical octupole Love number 

h 

0.0151 

Horizontal octupole Love number 

9 

980.665 

Surface acceleration due to gravity (cm/sec 2 ) 
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SECTION 8 


POSSIBLE IMPROVEMENTS TO THE CURRENT MODEL 


This section lists areas in which the current model can be improved. 

General Relativity: 

Variations of the Earth’s gravitational potential must be taken into account in defining 
proper lengths. This correction is estimated by Thomas (1991) to amount to 0.2 cm for a 
10,000 km baseline. 

Second-order effects have not been carefully investigated, and could possibly contribute at 
the picosecond level. 

Earth Orientation Models: 

There are short-period deficiencies in the present IAU models for the orientation of the 
Earth in space that may be as large as 1 to 2 milliarcseconds, and longer-term deficiencies 
of the order of 1 milliarcsecond per year (3 cm at one Earth radius). VLBI measurements 
made during the past few years indicate the need for revisions of this order of the annual 
nutation terms and the precession constant [Eubanks et al. (1985), Herring et al. (1986)]. 
The 18.6-year term in the IAU nutation series may also be in error, and present data spans 
are just approaching durations long enough to separate it from precession. To provide an 
improved nutation model, we have implemented a MODEST option to use the amplitudes 
of Zhu et al. (discussed in Section 2.6.2. 1). This will constitute a temporarily better model 
of the annual and semiannual nutations until the IAU series is officially revised. 

Tidal Effects: 

Ocean tides affect UTl, necessitating revisions and additional terms in the Yoder short- 
period UTl correction series (Brosche et al. } 1989). 

Antenna Deformation: 

Gravity loading and temperature variations may cause variations in the position of the 
reference point of a large antenna that are as large as 1 cm. Liewer (1986) presents evidence 
that these effects cause systematic errors and that their dependence on antenna orientation 
and ambient temperature may be modeled. 

Antenna Alignment: 

Hour angle misalignment of the order of 1 arc minute can cause 1 mm delay effects for DSN 
HA-Dec antennas with 7-m axis offsets. 

Subreflector Focusing: 

For DSN 70-m Cassegrain antennas, allowing the subreflector to slew in order to maintain 
focus changes the path delay by «8 cm over the 6° — 90° elevation range. Simulations 
(Jacobs, 1987) show that this effect is almost entirely absorbed by the clock epoch and local 
station vertical coordinate parameters. For baselines between two 70-m antennas, this causes 
a potential error of up to 12 cm in length. Presently, this effect can be modeled as a site 
vector relating fixed and slewed antenna positions; it may be more convenient to introduce 
a “slew flag” in the data to model it automatically. 

Phase Delay Rate: 

Rather than modeling the delay rates as finite differences of model delays, direct analytic 
expressions for derivatives of delays could be implemented. This would eliminate questions 
concerning the choice of the time difference A discussed in section 6. Care must be exercised, 
however, to ensure consistency between definitions of modeled and observed delay rates. 
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APPENDIX A 


NUTATION MODELS 


The three nutation series available in MODEST are tabulated here: Table A. I gives the standard 1980 
IAU series; Tables A.II, A.III, and A.IV contain the results of Zhu et al. (1990); for completeness, the 
old (Woolard) nutation series is given in Table A.V. 


Table A.I 


1980 IAU Theory of Nutation 


Index Period 
j (days) 


Argument coefficient 


fCyi Kj 2 *Cj3 


( 0 ". 0001 ) 


B 0 j Bij 
( 0 ". 0001 ) 


6798.4 

3399.2 

1305.5 

1095.2 
1615,7 
3232.9 

6786.3 

943.2 
182.6 

365.3 
121.7 
365.2 

177 8 

9 



-171996 

- 174.2 

92025 

2062 

0.2 

-895 

46 

0.0 

-24 

11 

0.0 

0 

-3 

0.0 

1 

-3 

0.0 

0 

-2 

0.0 

1 

1 

0.0 

0 

-13187 

- 1.6 

5736 

1426 

- 3.4 

54 

-517 

1.2 

224 

217 

- 0.5 

-95 

129 

0.1 

o 

i 

48 

0.0 

1 

-22 

0.0 

0 

17 

- 0.1 

0 

-15 

0.0 

9 

-16 

0.1 

7 

-12 

0.0 

6 

-6 

0.0 

3 

-5 

0.0 

3 


-2 0.0 

0 0.0 

0 0.0 

0 0.0 


0 0.0 














Table A. II 


Zhu et al. Theory of Nutation: 1980 IAU Terms 


Index 

j 

Period 

(days) 

kj i 

Argument coefficient 
kj 2 kj 3 & j '4 kj'5 

j4o/ Aij 

(0".00001) 

B oi Bij 
(0". 00001) 

1 

6798.38 

0 

0 

0 

0 

1 

-1720618 

-1743 

920530 

90 

2 

3399.19 

0 

0 

0 

0 

2 

20743 

2 

-8975 

5 

3 

1305.48 

-2 

0 

2 

0 

1 

460 

1 

-243 

0 

4 

1095.18 

2 

0 

-2 

0 

0 

110 

0 

1 

0 

5 

1615.75 

-2 

0 

2 

0 

2 

-31 

0 

14 

0 

6 

3232.86 

1 

-1 

0 

-1 

0 

-33 

0 

0 

0 

7 

6786.32 

0 

-2 

2 

-2 

1 

-15 

0 

8 

0 

8 

943.23 

2 

0 

-2 

0 

1 

7 

0 

-4 

0 

9 

182.62 

0 

0 

2 

-2 

2 

-131720 

-16 

57320 

-31 

10 

365.26 

0 

1 

0 

0 

0 

14735 

-35 

719 

-2 

11 

121.75 

0 

1 

2 

-2 

2 

-5176 

12 

2247 

-7 

12 

365.22 

0 

-1 

2 

-2 

2 

2161 

-5 

-961 

3 

13 

177.84 

0 

0 

2 

-2 

1 

1293 

1 

-699 

0 

14 

205.89 

2 

0 

0 

-2 

0 

479 

0 

5 

0 

15 

173.31 

0 

0 

2 

-2 

0 

-218 

0 

-1 

0 

16 

182.63 

0 

2 

0 

0 

0 

168 

-1 

2 

0 

17 

386.00 

0 

1 

0 

0 

1 

-140 

0 

86 

0 

18 

91.31 

0 

2 

2 

-2 

2 

-158 

1 

69 

0 

19 

346.64 

0 

-1 

0 

0 

1 

-127 

0 

64 

0 

20 

199.84 

-2 

0 

0 

2 

1 

-58 

0 

30 

0 

21 

346.60 

0 

-1 

2 

-2 

1 

-48 

0 

27 

0 

22 

212.32 

2 

0 

0 

-2 

1 

41 

0 

-22 

0 

23 

119.61 

0 

1 

2 

-2 

1 

36 

0 

-20 

0 

24 

411.78 

1 

0 

0 

-1 

0 

-43 

0 

-6 

0 

25 

131.67 

2 

1 

0 

-2 

0 

11 

0 

0 

0 

26 

169.00 

0 

0 

-2 

2 

1 

9 

0 

-4 

0 

27 

329.79 

0 

1 

-2 

2 

0 

-9 

0 

0 

0 

28 

409.23 

0 

1 

0 

0 

2 

7 

0 

-3 

0 

29 

388.27 

-1 

0 

0 

1 

1 

9 

0 

-4 

0 

30 

117.54 

0 

1 

2 

-2 

0 

-6 

0 

0 

0 

31 

13.66 

0 

0 

2 

0 

2 

-22824 

-2 

9806 

-5 

32 

27.55 

1 

0 

0 

0 

0 

7122 

1 

-70 

0 

33 

13.63 

0 

0 

2 

0 

1 

-3885 

-4 

2011 

0 

34 

9.13 

1 

0 

2 

0 

2 

-3023 

0 

1293 

-1 

35 

31.81 

1 

0 

0 

-2 

0 

-1572 

0 

-13 

0 

36 

27.09 

-1 

0 

2 

0 

2 

1238 

0 

-535 

0 

37 

14.77 

0 

0 

0 

2 

0 

635 

0 

-13 

0 

38 

27.67 

1 

0 

0 

0 

1 

633 

1 

-332 

0 

39 

27.44 

-1 

0 

0 

0 

1 

-580 

-1 

315 

0 

40 

9.56 

-1 

0 

2 

2 

2 

-598 

0 

256 

0 
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Table A. II cont. 


Zhu et al. Theory of Nutation: 1980 IAU Terms 


Index 

j 

Period 

(days) 

Argument coefficient 
kji kj 2 kjs 

fc /5 

(0".00001) 

Boj Bij 

(0".00001) 

41 

9.12 

1 

0 

2 

0 

1 

-517 

0 

265 

0 

42 

7.10 

0 

0 

2 

2 

2 

-386 

0 

165 

0 

43 

13.78 

2 

0 

0 

0 

0 

293 

0 

-6 

0 

44 

23.94 

1 

0 

2 

-2 

2 

286 

0 

-124 

0 

45 

6.86 

2 

0 

2 

0 

2 

-311 

0 

132 

0 

46 

13.61 

0 

0 

2 

0 

0 

259 

0 

-5 

0 

47 

26.98 

-1 

0 

2 

0 

1 

205 

0 

-107 

0 

48 

31.96 

-1 

0 

0 

2 

1 

152 

0 

-80 

0 

49 

31.66 

1 

0 

0 

-2 

1 

-129 

0 

70 

0 

50 

9.54 

-1 

0 

2 

2 

1 

-103 

0 

53 

0 

51 

34.85 

1 

1 

0 

-2 

0 

-74 

0 

-1 

0 

52 

13.17 

0 

1 

2 

0 

2 

76 

0 

-33 

0 

53 

14.19 

0 

-1 

2 

0 

2 

-71 

0 

31 

0 

54 

5.64 

1 

0 

2 

2 

2 

-77 

0 

32 

0 

55 

9.61 

1 

0 

0 

2 

0 

66 

0 

-3 

0 

56 

12.81 

2 

0 

2 

-2 

2 

65 

0 

-28 

0 

57 

14.80 

0 

0 

0 

2 

1 

-64 

0 

33 

0 

58 

7.09 

0 

0 

2 

2 

1 

-66 

0 

34 

0 

59 

23.86 

1 

0 

2 

-2 

1 

58 

0 

-30 

0 

60 

14.73 

0 

0 

0 

-2 

1 

-50 

0 

28 

0 

61 

29.80 

1 

-1 

0 

0 

0 

47 

0 

-1 

0 

62 

6.85 

2 

0 

2 

0 

1 

-53 

0 

27 

0 

63 

15.39 

0 

1 

0 

-2 

0 

-44 

0 

-1 

0 

64 

26.88 

1 

0 

-2 

0 

0 

41 

0 

1 

0 

65 

29.53 

0 

0 

0 

1 

0 

-40 

0 

1 

0 

66 

25.62 

1 

1 

0 

0 

0 

-34 

0 

1 

0 

67 

9.11 

1 

0 

2 

0 

0 

34 

0 

-1 

0 

68 

9.37 

1 

-1 

2 

0 

2 

-29 

0 

12 

0 

69 

9.81 

-1 

-1 

2 

2 

2 

-29 

0 

12 

0 

70 

13.75 

-2 

0 

0 

0 

1 

-23 

0 

13 

0 

71 

5.49 

3 

0 

2 

0 

2 

-29 

0 

12 

0 

72 

7.24 

0 

-1 

2 

2 

2 

-26 

0 

11 

0 

73 

8.91 

1 

1 

2 

0 

2 

25 

0 

-10 

0 

74 

32.61 

-1 

0 

2 

-2 

1 

-20 

0 

11 

0 

75 

13.81 

2 

0 

0 

0 

1 

22 

0 

-11 

0 

76 

27.78 

1 

0 

0 

0 

2 

-20 

0 

8 

0 

77 

9.18 

3 

0 

0 

0 

0 

16 

0 

-1 

0 

78 

9.34 

0 

0 

2 

1 

2 

16 

0 

-7 

0 

79 

27.33 

-1 

0 

0 

0 

2 

14 

0 

-6 

0 

80 

10.08 

1 

0 

0 

-4 

0 

-14 

0 

-1 

0 


71 



Table A. II cont. 


Zhu et al. Theory of Nutation: 1980 IAU Terms 


Index 

j 

Period 

(days) 

Argument coefficient 
kji 3 k/4 

k}5 

Aoj -4 ly 

(0".00001) 

Boj B l} 

(0".00001) 

81 

14.63 

-2 

0 

2 

2 

2 

13 

0 

-6 

0 

82 

5.80 

-1 

0 

2 

4 

2 

-15 

0 

6 

0 

83 

15.91 

2 

0 

0 

-4 

0 

-13 

0 

0 

0 

84 

22.47 

1 

1 

2 

-2 

2 

13 

0 

-5 

0 

85 

5.64 

1 

0 

2 

2 

1 

-13 

0 

7 

0 

86 

7.35 

-2 

0 

2 

4 

2 

-12 

0 

5 

0 

87 

9.06 

-1 

0 

4 

0 

2 

11 

0 

-5 

0 

88 

29.26 

1 

-1 

0 

-2 

0 

9 

0 

0 

0 

89 

12.79 

2 

0 

2 

-2 

1 

10 

0 

-5 

0 

90 

4.68 

2 

0 

2 

2 

2 

-11 

0 

5 

0 

91 

9.63 

1 

0 

0 

2 

1 

-10 

0 

5 

0 

92 

12.66 

0 

0 

4 

-2 

2 

9 

0 

-4 

0 

93 

8.75 

3 

0 

2 

-2 

2 

9 

0 

-4 

0 

94 

23.77 

1 

0 

2 

-2 

0 

-7 

0 

0 

0 

95 

13.41 

0 

1 

2 

0 

1 

8 

0 

-4 

0 

96 

35.03 

-1 

-1 

0 

2 

1 

7 

0 

-4 

0 

97 

13.58 

0 

0 

-2 

0 

1 

-6 

0 

3 

0 

98 

25.42 

0 

0 

2 

-1 

2 

-7 

0 

3 

0 

99 

14.19 

0 

1 

0 

2 

0 

-6 

0 

0 

0 

100 

9.53 

1 

0 

-2 

-2 

0 

-6 

0 

0 

0 

101 

14.16 

0 

-1 

2 

0 

1 

-7 

0 

3 

0 

102 

34.67 

1 

1 

0 

-2 

1 

-6 

0 

3 

0 

103 

32.76 

1 

0 

-2 

2 

0 

-6 

0 

0 

0 

104 

7.13 

2 

0 

0 

2 

0 

6 

0 

0 

0 

105 

4.79 

0 

0 

2 

4 

2 

-7 

0 

3 

0 

106 

27.32 

0 

1 

0 

1 

0 

5 

0 

0 

0 


Table A. Ill 


Zhu et al. Theory of Nutation: Out-of-Phase Terms 


Index 

j 

Period 

(days) 

fc/i 

Argument coefficient 

kj3 fcj'4 

kjs 

A^j B'2j 

(0".00001) 

1 

6798.38 

0 

0 

0 

0 

1 

221 112 

2 

182.62 

0 

0 

2 

-2 

2 

-153 -61 

3 

365.26 

0 

1 

0 

0 

0 

-55 22 

4 

13.66 

0 

0 

2 

0 

2 

-5 -2 
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Table A. IV cont. 


Zhu et al. Theory of Nutation: Planetary Terms 


Index 

j 

Period 

(days) 

Argument coefficient 
kji kj 2 kj 3 kj 4 


Ao} Boj 

(0".00001) 

41 

9.33 

0 

0 

2 

1 

1 

3 

-1 

42 

9.35 

1 

-1 

2 

0 

1 

-4 

2 

43 

9.60 

-1 

0 

0 

-2 

1 

-4 

3 

44 

10.07 

1 

0 

0 

-4 

1 

-1 

1 

45 

10.10 

-1 

0 

0 

4 

1 

-2 

1 

46 

10.37 

-1 

-1 

0 

4 

0 

1 

0 

47 

12.38 

2 

1 

2 

-2 

2 

3 

-1 

48 

12.64 

0 

0 

4 

-2 

1 

2 

-1 

49 

13.22 

1 

0 

2 

-1 

2 

-3 

1 

50 

13.28 

2 

1 

0 

0 

0 

-3 

0 

51 

13.63 

0 

0 

2 

0 

1 

-1 

0 

52 

13.69 

0 

0 

2 

0 

3 

2 

0 

53 

14.22 

0 

1 

0 

2 

1 

2 

-1 

54 

14.25 

1 

0 

0 

1 

0 

-3 

0 

55 

14.32 

2 

-1 

0 

0 

0 

4 

0 

56 

14.60 

-2 

0 

2 

2 

1 

2 

-1 

57 

14.70 

0 

0 

0 

-2 

2 

1 

-1 

58 

15.35 

0 

1 

0 

-2 

1 

-3 

2 

59 

15.42 

0 

-1 

0 

2 

1 

-2 

1 

60 

15.87 

2 

0 

0 

-4 

1 

-1 

1 

61 

15.94 

-2 

0 

0 

4 

1 

1 

-1 

62 

16.06 

0 

-2 

0 , 

2 

0 

2 

0 

63 

16.10 

0 

0 

2 

-4 

1 

-1 

1 

64 

22.40 

1 

1 

2 

-2 

1 

3 

-1 

65 

25.22 

-1 

1 

2 

0 

2 

4 

-2 

66 

25,53 

-1 

-1 

0 

0 

1 

2 

-1 

67 

25.72 

1 

1 

0 

0 

1 

-3 

2 

68 

26.77 

1 

0 

-2 

0 

1 

3 

-1 

69 

27.32 

0 

0 

1 

0 

1 

-2 

0 

70 

29.26 

-1 

-1 

2 

0 

2 

-2 

1 

71 

29.39 

-1 

1 

0 

2 

1 

-1 

1 

72 

29.40 

0 

0 

0 

-1 

1 

3 

-2 

73 

29.66 

0 

0 

0 

1 

1 

-4 

2 

74 

29.67 

-1 

1 

0 

0 

1 

-2 

2 

75 

31.52 

1 

0 

0 

-2 

2 

3 

-1 

76 

32.11 

-1 

0 

0 

2 

2 

-4 

2 

77 

32.45 

-1 

0 

2 

-2 

2 

3 

-1 

78 

35.80 

-1 

1 

2 

-2 

1 

-1 

1 

79 

38.52 

-1 

-2 

0 

2 

0 

3 

0 

80 

38.74 

1 

0 

2 

-4 

1 

-4 

2 
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Table A. IV cont. 


Zhu et al. Theory of Nutation: Planetary Terms 


Period 


Argument coefficient 


A 


(days) 


121.75 

129.17 
177.85 

219.17 
285.41 
297.91 
313.04 
329.82 
438.33 
471.95 
507.16 
552.62 

2266.13 

6159.14 
.74 
.86 
.58 
.73 
.82 
>.64 
>.73 
>.89 
>.95 
>.97 
>.98 
r .22 
f .50 
r .54 
5.94 
).10 
).20 
3.30 
3.37 
3.89 
3.08 
2.35 

2.71 
2.76 
3.49 

3.72 





Table A. IV cont. 


Zhu et al. Theory of Nutation: Planetary Terms 


Index 

j 

Period 

(days) 

k 3 i 

Argument coefficient 

kj 2 kj 3 kj 4 

kj 5 

121 

13.83 

2 

0 

0 

0 

2 

122 

14.13 

-1 

0 

2 

1 

2 

123 

14.16 

0 

1 

0 

2 

-1 

124 

14.76 

0 

-2 

2 

0 

2 

125 

14.93 

2 

0 

-2 

2 

-1 

126 

15.24 

-2 

-1 

2 

2 

2 

127 

15.31 

-1 


0 

3 

0 

128 


-2 

-1 

0 

4 

0 

129 

23.43 

-1 

0 

4 

-2 

2 

130 

23.94 

1 

2 

0 

0 

0 

131 

25.13 

-1 

1 

2 


-1 

132 

25.32 

0 

0 

2 

-1 

1 

133 

25.52 

1 

-1 

2 

-2 

1 

134 

25.62 

1 

-1 

2 

-2 

2 

135 

25.83 

2 

0 


-1 

0 

136 

27.09 

-1 

2 

0 

2 

0 

137 

27.32 

0 

-1 

2 

-1 

2 

138 

28.15 

3 


-2 

0 

-1 

139 

29.14 

-1 

-1 

2 

0 

1 

140 

29.14 

-1 

1 

0 

2 

-1 

141 

31.06 

-3 

0 


2 

1 

142 

32.45 

1 

-2 

0 

0 

0 

143 

34.48 

-2 

0 

0 

3 

0 

144 

37.62 

-3 

0 

0 

4 

0 

145 

38.52 

-1 

0 

-2 

4 

-2 

146 

38.96 

-1 

0 

-2 

4 

0 

147 

43.06 

-1 

-1 

-2 

4 

-2 

148 

43.34 

1 

1 

2 

-4 

1 

149 

90.10 

0 


2 

-2 

1 

150 

96.78 

2 

0 

2 

-4 

2 

151 

134.27 

2 

1 

0 

-2 

1 

152 

156.52 

-2 


4 

-2 

2 

153 

164.08 

-2 

2 

2 

0 

2 

154 

187.67 

0 

2 

0 

0 

1 

155 

193.56 


-1 

2 

-3 

2 

156 

235.96 


0 

2 

2 

2 


Aoy B$j 
( 0 ". 00001 ) 


-1 0 

1 0 

-1 0 

-1 0 

1 0 

1 

-1 
1 

-1 
-1 
1 

-1 
-1 
-1 
1 

-1 
1 
1 

-1 
-1 
1 

1 0 

-1 0 

1 0 

-1 0 

-1 0 

1 0 

-1 0 

1 0 

-1 0 

1 0 

-1 0 

-1 0 

-1 0 

1 0 

-1 0 
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Table A.V 


Woolard Theory of Nutation 


Index 

j 

Period 

( days ) 

Argument coefficient 
kji kj 2 kj 3 i /4 

fc/5 

•^Oj A lj 

(o".oooi) 5 

Boj Bij 
( 0 ". 0001 ) 

1 

6798.4 

0 

0 

0 

0 

1 

-172327 - 

173.7 

92100 

9.1 

2 

3399.2 

0 

0 

0 

0 

2 

2088 

0.2 

-904 

0.4 

3 

1305.5 

-2 

0 

2 

0 

1 

45 

0.0 

-24 

0.0 

4 

1095.2 

2 

0 

-2 

0 

0 

10 

0.0 

0 

0.0 

5 

1615.7 

-2 

0 

2 

0 

2 

-3 

0.0 

2 

0.0 

6 

3232.9 

1 

-1 

0 

-1 

0 

-2 

0.0 

0 

0.0 

7 

6786.3 

0 

-2 

2 

-2 

1 

-4 

0.0 

2 

0.0 

8 

182.6 

0 

0 

2 

-2 

2 

-12729 

- 1.3 

5522 

- 2.9 

9 

365.3 

0 

1 

0 

0 

0 

1261 

- 3.1 

0 

0.0 

10 

121.7 

0 

1 

2 

-2 

2 

-497 

1.2 

216 

- 0.6 

11 

365.2 

0 

-1 

2 

-2 

2 

214 

- 0.5 

-93 

0.3 

12 

177.8 

0 

0 

2 

-2 

1 

124 

0.1 

-66 

0.0 

13 

205.9 

2 

0 

0 

-2 

0 

45 

0.0 

0 

0.0 

14 

173.3 

0 

0 

2 

-2 

0 

-21 

0.0 

0 

0.0 

15 

182.6 

0 

2 

0 

0 

0 

16 

- 0.1 

0 

0.0 

16 

386.0 

0 

1 

0 

0 

1 

-15 

0.0 

8 

0.0 

17 

91.3 

0 

2 

2 

-2 

2 

-15 

0.1 

7 

0.0 

18 

346.6 

0 

-1 

0 

0 

1 

-10 

0.0 

5 

0.0 

19 

199.8 

-2 

0 

0 

2 

1 

-5 

0.0 

3 

0.0 

20 

346.6 

0 

-1 

2 

-2 

1 

-5 

0.0 

3 

0.0 

21 

212.3 

2 

0 

0 

-2 

1 

4 

0.0 

-2 

0.0 

22 

119.6 

0 

1 

2 

-2 

1 

3 

0.0 

-2 

0.0 

23 

411.8 

1 

0 

0 

-1 

0 

-3 

0.0 

0 

0.0 

24 

13.7 

0 

0 

2 

0 

2 

-2037 

- 0.2 

884 

- 0.5 

25 

27.6 

1 

0 

0 

0 

0 

675 

0.1 

0 

0.0 

26 

13.6 

0 

0 

2 

0 

1 

-342 

- 0.4 

183 

0.0 

27 

9.1 

1 

0 

2 

0 

2 

-261 

0.0 

113 

- 0.1 

28 

31.8 

1 

0 

0 

-2 

0 

-149 

0.0 

0 

0.0 

29 

27.1 

-1 

0 

2 

0 

2 

114 

0.0 

-50 

0.0 

30 

14.8 

0 

0 

0 

2 

0 

60 

0.0 

0 

0.0 

31 

27.7 

1 

0 

0 

0 

1 

58 

0.0 

-31 

0.0 

32 

27.4 

-1 

0 

0 

0 

1 

-57 

0.0 

30 

0.0 

33 

9.6 

-1 

0 

2 

2 

2 

-52 

0.0 

22 

0.0 

34 

9.1 

1 

0 

2 

0 

1 

-44 

0.0 

23 

0.0 

35 

7.1 

0 

0 

2 

2 

2 

-32 

0.0 

14 

0.0 
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Table A.V confc. 


Woolard Theory of Nutation 


Index 

j 

Period 

(days) 

Argument coefficient 
kji kj 2 Ary 3 fcy 4 

Arys 

-^oy ^iy 

(0".0001) 

Soy Siy 
(0".0001) 

36 

13.8 

2 

0 

0 

0 

0 

28 

0.0 

0 

0.0 

37 

23.9 

1 

0 

2 

-2 

2 

26 

0.0 

-11 

0.0 

38 

6.9 

2 

0 

2 

0 

2 

-26 

0.0 

11 

0.0 

39 

13.6 

0 

0 

2 

0 

0 

25 

0.0 

0 

0.0 

40 

27.0 

-1 

0 

2 

0 

1 

19 

0.0 

-10 

0.0 

41 

32.0 

-1 

0 

0 

2 

1 

14 

0.0 

-7 

0.0 

42 

31.7 

1 

0 

0 

-2 

1 

-13 

0.0 

7 

0.0 

43 

9.5 

-1 

0 

2 

2 

1 

-9 

0.0 

5 

0.0 

44 

34.8 

1 

1 

0 

-2 

0 

-7 

0.0 

0 

0.0 

45 

13.2 

0 

1 

2 

0 

2 

7 

0.0 

-3 

0.0 

46 

14.2 

0 

-1 

2 

0 

2 

-6 

0.0 

3 

0.0 

47 

5.6 

1 

0 

2 

2 

2 

-6 

0.0 

3 

0.0 

48 

9.6 

1 

0 

0 

2 

0 

6 

0.0 

0 

0.0 

49 

12,8 

2 

0 

2 

-2 

2 

6 

0.0 

-2 

0.0 

50 

14.8 

0 

0 

0 

2 

1 

-6 

0.0 

3 

0.0 

51 

7.1 

0 

0 

2 

2 

1 

-5 

0.0 

3 

0.0 

52 

23.9 

1 

0 

2 

-2 

1 

5 

0.0 

-3 

0.0 

53 

14.7 

0 

0 

0 

-2 

1 

-5 

0.0 

3 

0.0 

54 

29.8 

1 

-1 

0 

0 

0 

4 

0.0 

0 

0.0 

55 

6.9 

2 

0 

2 

0 

1 

-4 

0.0 

2 

0.0 

56 

15.4 

0 

1 

0 

-2 

0 

-4 

0.0 

0 

0.0 

57 

26.9 

1 

0 

-2 

0 

0 

4 

0.0 

0 

0.0 

58 

29.5 

0 

0 

0 

1 

0 

-4 

0.0 

0 

0.0 

59 

25.6 

1 

1 

0 

0 

0 

-3 

0.0 

0 

0.0 

60 

9.1 

1 

0 

2 

0 

0 

3 

0.0 

0 

0.0 

61 

9.4 

1 

-1 

2 

0 

2 

-3 

0.0 

0 

0.0 

62 

9.8 

-1 

-1 

2 

2 

2 

-2 

0.0 

0 

0.0 

63 

13.7 

-2 

0 

0 

0 

1 

-2 

0.0 

0 

0.0 

64 

5.5 

3 

0 

2 

0 

2 

-2 

0.0 

0 

0,0 

65 

7.2 

0 

-1 

2 

2 

2 

-2 

0.0 

0 

0.0 

66 

8.9 

1 

1 

2 

0 

2 

2 

0.0 

0 

0.0 

67 

32.6 

-1 

0 

2 

-2 

1 

-2 

0.0 

0 

0.0 

68 

13.8 

2 

0 

0 

0 

1 

2 

0.0 

0 

0.0 

69 

27.8 

1 

0 

0 

0 

2 

-2 

0.0 

0 

0.0 
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APPENDIX B 


GLOSSARY OF “MODEST” PARAMETERS 

For the convenience of users of MODEST, Table B.I identifies the names of adjustable parameters in 
the code with the notation of this document. Brief definitions and either references to equations (in 
parentheses) or sections (no parentheses) are also given. 

Table B.I 

Glossary of MODEST Parameters 


Parameter 

MODEST name 

Definition 

Reference 

r ap 

RSPINAX aaaaaaaa 

Cylindrical 

(2.38) 

A 

LQNGTUD aaaaaaaa 

station 

(2.39) 

z 

POLPROJ aaaaaaaa 

coordinates 

(2.40) 


DRSP/DT aaaaaaaa 

Time rates of 

(2.38) 

A 

DLON/DT aaaaaaaa 

change of 

(2.39) 

z 

DPOL/DT aaaaaaaa 

stn. coords. 

(2.40) 

X 

X aaaaaaaa 

Cartesian 

(2.41) 

y 

Y aaaaaaaa 

station 

(2.42) 

z 

Z aaaaaaaa 

coordinates 

(2.43) 

X 

DX/DT aaaaaaaa 

Time rates of 

(2.41) 

y 

DY/DT aaaaaaaa 

change of 

(2.42) 

z 

DZ/DT aaaaaaaa 

stn. coords. 

(2.43) 

l 

AXISOFF aaaaaaaa 

Antenna offset 

(2.172) 

h, l 

♦LOVE # aaaaaaaa 

Love numbers 

(2.51) to (2.53) 


TIDEPHZ aaaaaaaa 

Tide lag 

(2.48) 

~1pPN 

GEN REL GAMMA FACTOR 

PPN gamma 

(2.16) 

a 

RIGHT ASCEN . ssssssssssss 

Source RA 

(2.199) 

6 

DECLINATION ssssssssssss 

Source dec. 

(2.199) 

a 

DRASCEN/DT ssssssssssss 

Time rates of 

(2.85) 

6 

DDECLIN/DT ssssssssssss 

change of RA, dec. 

(2.86) 

©1,2 

t POLE MOTION 

Pole position 

(2.90), (2.91) 

UT\ - UTC 

UT1 MINUS UTC 

UTl - UTC 

2.6.1 


aaaaaaaa station name 

ssssssssssss source name 
* V or H 
t X or Y 
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Table B.I cont. 

Glossary of MODEST Parameters 


Parameter 

MODEST name 

Definition 

Reference 


t AXIS TWEAK OFFSET 

Perturbation 

(2.138) 


X AXIS TWEAK RATE 

coefficients 

(2.138) 

Pls 

LUNI-SOLAR PRECESSION 

Precession 

(2.128) 

PPL 

PLANETARY PRECESSION 

constants 

(2.128) 

Aoj 

NUTATION AMPLTD PSI cjjj 

Nutation 

(2.113) to 

*xi 

A Wj 
B 'oy 

Bxs 

NUTATION AMPLTD PSITcjjj 
NUTATION AMPLTD PSIA 
NUTATION AMPLTD EPS cjjj 
NUTATION AMPLTD EPSTcjjj 
NUTATION AMPLTD EPSA 

amplitudes 

(2.118) 

Tel 

C EPOCH aaaaaaaa 

Coefficients 

(3.1) 

T C 2 

C RATE aaaaaaaa 

in clock 

(3.1) 

T c 3 

DCRAT/DTaaaaaaaa 

model for 

(3.1) 

^o4 

F OFFSETaaaaaaaa 

delay and 

(3.2) 

Teh 

F DRIFT aaaaaaaa 

delay rate 

(3.2) 

PZdry 

DRYZTROPaaaaaaaa 

Dry zenith delay 

(4.3) 

PZv'f 

WETZTROPaaaaaaaa 

Wet zenith delay 

(4.3) 

PZdry 

DDTRP/DTaaaaaaaa 

Zenith delay 

(4.4) 

PZuet 

DWTRP/DTaaaaaaaa 

time rates 

(4.4) 

Adry 

DRYZMAPAaaaaaaaa 

Chao map 

(4.7) to 

&dry 

DRYZMAPBaaaaaaaa 

parameters 

(4.11) 

P 

DRYMAPSGaaaaaaaa 

Lanyi map 
parameter 

(4.30) 

To 

SURFTEMPaaaaaaaa 

CfA map surface 
temperature 

(4.36) 

le add 

Z TECADDaaaaaaaa 

Zenith electron 
content 

(5.23) 


t X, Y, or Z 
c 

JJJ 

aaaaaaaa 


component: S, C for sine, cosine 
1980 IAU series term number 
station name 
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